{
 "metadata": {
  "name": "Convergent-divergent nozzle"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Convergent-divergent nozzle"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Statement**: Determine the steady flow within the convergent-divergent nozzle, geometrically defined by\n",
      "\n",
      "$$A(x) = \\begin{cases}\n",
      "2 x^3 - 3 x^2 + 2 & 0 \\leq x \\le 1 \\\\\\\\\n",
      "-\\frac{3}{8} x^3 + \\frac{9}{4} x^2 - \\frac{27}{8} x + \\frac{5}{2} & 1 \\leq x \\le 3\n",
      "\\end{cases}$$\n",
      "\n",
      "for the operating condition in which the nozzle discharges in the atmosphere and is supplied with a stagnation pressure P_0 = 2 atm and a stagnation temperature T0 = 300 \u00b0C."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Algorithm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Temptative algorithm: given nozzle geometry $A(x)$, ambient pressure $p_A$ and inlet stagnation conditions $p_0$ and $T_0$,\n",
      "\n",
      "1. We are discharging to the atmosphere, therefore the back pressure $p_B$ is the given ambient pressure.\n",
      "2. Supose chocked flow, hence minimal area $A_t$ is critical area $A^*$.\n",
      "3. Compute $M_e$ both subsonic and supersonic, using the Mach number-area ratio relation.\n",
      "3. Compute subsonic $p_{e3}$ and supersonic $p_{e6}$ pressure at the exit using isentropic Mach number relations and inlet stagnation pressure.\n",
      "4. Compare $p_B$ with $p_{e3}$ and $p_{e6}$:\n",
      "    5. If $p_{e3} \\lt p_B$, fully subsonic flow, no sonic conditions at the throat.\n",
      "    6. If $p_B \\leq p_{e3}$, chocked flow.\n",
      "7. ..."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Computations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "\n",
      "# Import necessary packages\n",
      "import numpy as np\n",
      "import scipy as sp\n",
      "\n",
      "from scipy import constants, optimize\n",
      "\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n",
        "For more information, type 'help(pylab)'.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from skaero.gasdynamics import isentropic, shocks"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 0. Given conditions\n",
      "# Taken from exercises 2 and 18\n",
      "\n",
      "# Parameters\n",
      "gamma = g = 1.4\n",
      "R = 287.04  # J\u2009/ (kg \u00b7 K)\n",
      "p_amb = 1 * sp.constants.atm  # Pa\n",
      "\n",
      "# Inlet conditions\n",
      "p_0 = 2.0 * sp.constants.atm  # Pa\n",
      "T_0 = sp.constants.C2K(300.0)  # K\n",
      "\n",
      "# Longitudinal axis of the nozzle\n",
      "x = np.linspace(0, 3, 151)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Geometry of the nozzle\n",
      "def A(x):\n",
      "    \"\"\"Cross sectional area of the nozzle, in m^2.\n",
      "\n",
      "    Notice that this function also contains a couple of workarounds\n",
      "    to circunvent some odd behaviour of numpy.piecewise.\n",
      "    \"\"\"\n",
      "    def A_1(x):\n",
      "        return (2.0 * x ** 3 - 3.0 * x ** 2 + 2.0) / 100\n",
      "\n",
      "    def A_2(x):\n",
      "        return (-3.0 * x ** 3 / 8.0 + 9.0 * x ** 2 / 4.0 - 27.0 * x / 8.0 + 5.0 / 2.0) / 100\n",
      "\n",
      "    x = np.asarray(x, dtype=float)\n",
      "    # For avoiding np.piecewise bug\n",
      "    if not x.shape:\n",
      "        x = np.asarray([x])\n",
      "\n",
      "    result = np.piecewise(x, [(0.0 <= x) & (x < 1.0), (1.0 <= x) & (x <= 3.0)], [A_1, A_2])\n",
      "    if len(result) == 1:\n",
      "        result = result[0]\n",
      "\n",
      "    return result\n",
      "\n",
      "radius = np.sqrt(A(x) / np.pi) * 10  # dm\n",
      "\n",
      "# Plot nozzle shape\n",
      "nozzle = plt.fill_between(x, radius, -radius, facecolor=\"#cccccc\")\n",
      "plt.xlim(0, 3)\n",
      "plt.title(\"Nozzle\")\n",
      "plt.xlabel(\"x (m)\")\n",
      "plt.ylabel(\"Radius (dm)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.text.Text at 0x7fe7d80ade10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEXCAYAAABRWhj0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xlwk2d+B/CvbPm+bYwPSfi+ZFuSD3DILtcS4hgS1xDC\nEdgyObY0bTaz3W7atJNmm5bNJNvsNCGZzSTTaQjphqTZnVyQUEzwbfkWGOMbbOPb+L7kS3r7RxoV\ngx0LHzrs72fmHSO9j6Tfwwv+6n3e93lfkSAIAoiIiIxgY+4CiIjIejA0iIjIaAwNIiIyGkODiIiM\nxtAgIiKjMTSIiMhoDA0iM8rOzoZMJjN3GURGY2gQ3SE4OBh+fn4YHx83PPcf//Ef2LFjhxmrIrIM\nDA2iOej1erz55pvmLoPI4jA0iO4gEonwq1/9Cq+//jqGhobuWl9YWIiNGzfC09MTmzZtglqtBgCo\n1Wq4ubkZFkdHR4SEhAAAPD09Dc+7urrCxsYGN2/evOu9Ozo68Oijj2L9+vUIDQ3FW2+9tbKdJbpH\nDA2iOSQnJ2P79u14/fXXZz0/MDCAPXv24Be/+AX6+/vxy1/+Env27EF/fz82b96MkZERjIyMYGBg\nAPfddx8ef/xxAMDg4KBh3XPPPYetW7dCIpHMem+9Xo9HHnkECQkJ6OjowLfffos33ngDFy5cMFm/\niRbC0CCag0gkwr/8y7/grbfeQm9vr+H5c+fOITIyEkeOHIGNjQ0OHTqE6OhofPXVV7Ne//Of/xzu\n7u74zW9+M+v5Tz75BGfOnMGf/vQn2NrazlpXWlqK3t5evPjiixCLxQgJCcHTTz+Njz/+eOU6SnSP\nxOYugMhSxcbG4uGHH8arr76KmJgYAN8NHwUFBc1qFxQUhPb2dsPjd999F7m5uSguLp7VTqPR4Oc/\n/zkyMzPh4+Nz1+e1tLSgo6MDXl5ehud0Oh22bt26nN0iWhKGBtEPePnll5GYmIi//du/BQAEBgai\npaVlVpuWlhakpaUBAPLy8vDSSy+hoKAArq6uhjY9PT3Yu3cvfv/730OpVM75WTKZDCEhIaivr1+h\n3hAtHYeniH5AWFgYDh48iDfffBMikQi7d+9GfX09zpw5g5mZGXzyySeora3Fww8/jNbWVhw4cAAf\nfvghwsPDDe8xMzOD/fv34+jRo9i/f/+8n7Vp0ya4ubnht7/9LbRaLXQ6HaqqqlBWVmaKrhIZhaFB\ntICXXnrJMGfD29sbZ8+exe9+9zusW7cOr7/+Os6ePQtvb298++236OnpwaOPPmo4Uyo+Ph7t7e3I\nz8/HG2+8YXje3d0dbW1tAL47fgIAtra2OHv2LC5fvozQ0FD4+vriL/7iLzA8PGy2vhPdScSbMBER\nkbHMtqfx5JNPws/PD/Hx8fO2ee655xAREQGlUgmNRmPC6oiIaC5mC40nnngC58+fn3f9119/jcbG\nRjQ0NOC9997DM888Y8LqiIhoLmYLjS1btsw6tfBOX375JY4dOwYASElJweDgILq7u01VHhERzcFi\nT7ltb2+fdfVPqVSKtrY2+Pn5zWr3/UFEIiK6N4s5pG3RZ0/d2aH5AkIQhFW7/PrXvzZ7Dewb+8f+\nrb5lsSw2NCQSCVpbWw2P29ra7rpWDxERmZbFhkZ6ejpOnz4NACgqKoKnp+ddQ1NERGRaZjumcfjw\nYeTk5KC3txcymQwvv/wypqenAQDHjx/H7t278fXXXyM8PBwuLi54//33zVWqWW3fvt3cJayY1dw3\ngP2zdqu9f4tl9ZP7RCLRksbniIjWosX+7rTY4SkiIrI8DA0iIjIaQ4OIiIzG0CAiIqMxNIiIyGgM\nDSIiMprFXnuKiIiA8fFx9Pf3Y2BgAIODgxgcHMTQ0BCGh4cxPDyM0dFRDA8PY2xsDGNjY9BqtdBq\ntZicnMTk5CSmpqYwPT2NmZkZ6HS6JV9GhPM0iIhMSBAE9Pf3o6OjAx0dHejq6kJHRwd6enrQ1dWF\nnp4e9PT0oK+vDwMDAxAEAZ6ennB3d4e7uzvc3Nzg4uICZ2dnODs7w8nJybA4OjrCwcEBjo6OsLOz\ng729PcRisWGxsbGBjc13A0z79+9f1O9OhgYR0TIaHR1FS0vLrKW5uRltbW1ob29HZ2cnHBwcsH79\nevj6+sLHxwdeXl6Gn97e3vDy8oKnpyc8PT3h6Oi4InUmJyczNIiIVpogCOjs7ERjY6NhaWhowI0b\nN9DS0oKxsTFIJBIEBATAz88P69evh7+/P9avX294vFJBcC8YGkREy2hkZAR1dXWora01LPX19bh+\n/TqcnJwQFBQEiUSCwMBASCQSSKVSBAYGwsfHxyru88PQICJahMHBQVy7ds2wVFVVoba2Fv39/QgJ\nCUFwcDCkUimCgoKwYcMGbNiwAa6uruYue8kYGkREP2Bqagq1tbWorKzElStXcOXKFVy7dg2Dg4MI\nCwtDaGgogoODDT8DAgIMB41XI4YGEdH/GRgYwJUrV6DRaAxLY2MjJBIJwsPDERoaivDwcISFhSEw\nMHBVh8N8FhsanKdBRFbt1q1bKC8vR1lZGcrKynD58mX09vYiKioKERERiIyMxK5duxAaGmoRB6Ct\nHfc0iMhqDAwMoKysDCUlJSgpKUF5eTlGRkYgl8sRGRmJqKgoREVFYcOGDWty7+FerOnhqY8//hhi\nsRj29vZwcnKCi4sL3N3d4eHhgXXr1vHbBZEVmpiYwOXLl1FUVITi4mKUlJSgu7sbcrkc0dHRiI6O\nhlwuh1QqtYqzlSzNmg6NlJQUuLq6YmZmBlNTU9BqtYap9YODgxCLxfDz84O/vz9kMhmCg4MREhJi\nGNcMCgqCra2tubtCtGYJgoCmpiYUFRWhsLAQarUaNTU1CAkJQWxsLKKjoxEXF4fg4GD+X10mazo0\nTp06hbi4uDnXC4KAsbEx9PX1obe3F93d3ejq6jJM3b958yb6+/sRGhqK2NhYKBQKqFQqqFQqSCQS\nfoMhWgHj4+MoLS1FYWEhCgoKUFJSAhsbG8THxyM2NhZxcXGIiYnhKMEK4oHweYhEIri6usLV1RVB\nQUFzttFqtWhpacH169fR0NCA//mf/0FtbS0AIDExESkpKdi8eTNSUlLg7e1tyvKJVoWOjg4UFBQg\nPz8f+fn5qKmpQUREBOLi4vDjH/8YzzzzDPz9/c1dJhlh1e9pLJYgCLh16xaqq6tx7do1w0+JRIIt\nW7Zg27Zt2LZtG2Qy2bJ+LpG10+v1qK6uRl5eHvLy8lBQUICRkREolUrEx8dDoVBwL8ICcHhqmUNj\nLjMzM2hsbIRGo8GVK1dQUVEBDw8P/OQnP8GDDz6InTt3wsfHZ8XrILIkk5OTKC8vR25uLnJyclBU\nVAQPDw9DSCQkJCAoKIhDvRaGoWGC0LiTXq/H9evXUVpairKyMlRUVCAiIgK7d+/G7t27kZKSwoN2\ntOqMjo5CrVYjNzcX2dnZqKioQHBwMJRKJVQqFZRKJdatW2fuMmkBDA0zhMadpqenUVlZCbVajaKi\nIvT09CA1NRXp6el46KGH4Onpae4Sie7ZwMAA8vPzkZOTg+zsbNTU1CAmJsYQEgqFYlVci2mtYWhY\nQGjcqaurCwUFBSgsLER5eTmSk5Oxb98+ZGRk8FgIWayenh7DXkROTg6am5uhUCigUCiQmJiI2NhY\nODg4mLtMWiKGhgWGxu20Wi2Ki4uRm5uL3NxchIaGYv/+/di/fz8iIiLMXR6tYZ2dncjJyUFWVhZy\ncnLQ2dlpGGZKSkpCdHQ0xOJVf6LlmsPQsPDQuN3MzAwqKiqQlZWF7Oxs+Pn54cCBAzhw4ACioqLM\nXR6tcq2trbNCoq+vD4mJiVAqlUhMTERkZCSPxa0BDA0rCo3b6XQ6XLlyBZcuXcKlS5fg6+uLQ4cO\n4dChQwgPDzd3ebQKNDc3Izs72xASIyMjSExMREJCAhITExEWFsbrNK1BDA0rDY3b6fV6XL58GRcv\nXsSlS5cglUpx5MgRHDx4kMdAyCiCIOD69euz9iQmJiaQlJQElUqFxMREhIaG8vRXYmishtC43czM\nDMrLy5GZmYns7GzI5XIcOXIEBw4c4OmMZCAIAmpqagxnNuXk5AD47koGKpUKSUlJnCNBc2JorLLQ\nuN3U1BTUajUuXryI/Px8bN68GUePHsXevXt5quMao9PpUFlZaQiJgoICODk5ISEhwXDgmtdMI2Mw\nNFZxaNxufHwcubm5uHDhAjQaDR566CEcPXoUqampPA1yFZqcnERZWZnhFNiioiL4+vpCqVQiISEB\nCQkJvGYTLQpDY42Exu0GBwfx7bffIjMzE9evX8fevXtx5MgRbNu2jWe/WKnBwUGo1Wrk5eUhNzcX\nGo0GISEhhol0KpWKF82kZcHQWIOhcbuuri5cuHABFy9eRF9fHw4cOIAjR45g06ZNHKqwUIIgoKWl\nBQUFBcjLy0N+fj6am5sNl+j//tpNHIKklcDQWOOhcbvm5mZkZmYiMzMTMzMzOHjwIA4dOoSEhAQG\niBlNTEygoqLCcA8JtVoNvV5vCAelUsmJdGQyDA2Gxl0EQUB9fT0yMzPx7bffQiwW4+DBg3jssccY\nICtMEATcuHEDxcXFUKvVUKvVqK6uRmhoKOLj4xEXFweFQoGAgABuBzILhgZD4wcJgoDa2lpcvHgR\nWVlZEIlEePTRR7Fv3z5s3ryZk7uWqLOzE2VlZSgpKUFxcTHKy8vh6OgIuVwOuVyOuLg4yOVy3kOC\nLAZDg6FhtO/3QLKzs5Gbm4v+/n48/PDDyMjIwM6dO+Hi4mLuEi2WIAhoa2uDRqNBRUUFSktLUVFR\ngYmJCcO9rOVyOWJjYzmfhiwaQ4OhsWhtbW3Izc1FQUEBrl27hs2bN2PPnj3YvXs3wsPD1+zwiVar\nRW1tLSorK3H58mVoNBpcvXoVNjY2iImJQXh4OKKjoxEdHY3AwMA1+/dE1omhwdBYFqOjoygqKkJR\nURHUajXs7e3xwAMPYNeuXdi2bRsCAwPNXeKyGx8fR11dHWpqalBdXY2rV6+iuroabW1tkMlkiIiI\nQGhoKCIjIxEZGQkfHx8GBFk9hgZDY9kJgoDm5mYUFxdDo9GgvLwc3t7e+NGPfoQtW7Zg06ZNiI2N\ntYqzfUZHR9HU1IQbN26gsbER9fX1aGhoQH19Pfr6+rBhwwYEBwdjw4YNCA0NRWhoKIKCgmBnZ2fu\n0olWBEODobHi9Ho9bty4gStXruDq1auoqalBZ2cn5HK54TIWcrkc0dHR8Pf3N8m3cUEQMDQ0hO7u\nbnR0dKCzsxNtbW1oa2tDS0sLbt68idbWVoyNjUEikUAqlSIwMBBSqRQymQxBQUHw9/fnZEhacxga\nDA2zGB0dNXxjb2pqQnNzM5qamjA+Po6goCDIZDJIpVL4+flh/fr18PDwgJubG5ycnGBnZwdbW1vo\n9Xro9XpMT09jenoaU1NTmJiYwPj4OMbGxjA6Oorh4WEMDg5iYGAAAwMD6OvrQ19fH/r7+2Fvbw8f\nHx/4+voafvr6+sLf3x9+fn4ICAiAt7c3h5SIbrPY0LD8cQWyaK6uroZrIN1udHQUHR0duHXrFrq7\nuzEwMICWlhaMj49jfHwck5OT0Ol00Ov1EIlEEIlEEIvFEIvFsLOzg729PRwcHODo6AhHR0e4urpC\nIpEgJiYGbm5u8PDwgIeHBzw9PXkaK5EJMTRoRbi6uhoOHBPR6sEZXUREZDSGBhERGc2soXH+/HlE\nR0cjIiICr7322l3rs7Oz4eHhYRgzP3HihBmqJCKi75ntmIZOp8Ozzz6LixcvQiKRYOPGjUhPT0dM\nTMysdtu2bcOXX35ppiqJiOh2ZguNkpIShIeHIzg4GABw6NAhfPHFF3eFhjGnhH3++ecoLCwEACQl\nJSEpKWnZ6yUismbl5eUoLy9f8vuYLTTa29shk8kMj6VSKYqLi2e1EYlEKCwshFKphEQiweuvvw65\nXH7Xe2VkZHCeBhHRD7jzC/V77723qPcxW2gYM9EqMTERra2tcHZ2xjfffIOMjAzU19eboDoiIpqL\n2Q6ESyQStLa2Gh63trZCKpXOauPm5gZnZ2cAQFpaGqanp9Hf32/SOomI6P+ZLTSSk5PR0NCA5uZm\nTE1N4ZNPPkF6evqsNt3d3YZjGiUlJRAEAd7e3uYol4iIYMbhKbFYjLfffhupqanQ6XR46qmnEBMT\ng3fffRcAcPz4cfzxj3/EO++8A7FYDGdnZ3z88cfmKpeIiMALFhIRrUmLvWAhZ4QTEZHRGBpERGQ0\nhgYRERmNoUFEREZjaBARkdEYGkREZDSGBhERGY2hQURERmNoEBGR0RgaRERkNIYGEREZjaFBRERG\nY2gQEZHRGBpERGQ0hgYRERmNoUFEREZjaBARkdEYGkREZDSGBhERGY2hQURERmNoEBGR0cQ/tHJ6\nehoXLlxAbm4umpubIRKJEBQUhK1btyI1NRVi8Q++nIiIVpl59zT+9V//FRs3bsTZs2cRHR2NJ598\nEseOHUNUVBS++uorJCcn48SJE6aslYiIzGzeXQWlUokXX3wRIpHornVPPvkk9Ho9zp49u6LFERGR\nZZk3NNLT03/whTY2Ngu2ISKi1WXBgxKlpaV45ZVX0NzcjJmZGQCASCRCZWXlihdHRESWZcHQOHLk\nCF5//XXExcXBxoYnWxERrWULhoavry+HoYiICIARofHrX/8aTz31FB544AHY29sD+G54at++fSte\nHBERWZYFQ+ODDz5AXV0dZmZmZg1PMTSIiNaeBUOjrKwMtbW1c556S0REa8uCR7bvv/9+VFdXm6IW\nIiKycAvuaajVaqhUKoSEhMDBwQEAT7klIlqrFgyN8+fPm6IOIiKyAvOGRn9/PwDA3d3dZMUQEZFl\nmzc0EhMTIRKJIAgCbt68CS8vLwDAwMAAgoKC0NTUZLIiiYjIMsx7ILy5uRlNTU3YtWsXzp49i76+\nPvT19eHcuXPYtWuXKWskIiILseDZU2q1Grt37zY8TktLQ2Fh4YoWRURElmnBA+GBgYE4ceIEjh49\nCkEQ8NFHH0EikZiiNiIisjAL7mmcOXMGPT092Lt3L/bt24eenh6cOXPGFLUREZGFWXBPw8fHBydP\nnjRFLUREZOHm3dN48sknUVpaOu8Li4uL8cQTT6xIUUREZJnm3dP4m7/5G/zbv/0bioqKEBUVhYCA\nAAiCgK6uLtTV1eH+++/Hr371K1PWSkREZjZvaMTHx+P06dOYnJyERqNBS0sLRCIRgoKCoFQq4ejo\naMo6iYjIAix4TMPBwQH33Xcf7rvvPlPUQ0REFoz3byUiIqOZNTTOnz+P6OhoRERE4LXXXpuzzXPP\nPYeIiAgolUpoNBoTV0hERLe7p9DQ6XQYHh5elg/W6XR49tlncf78eVRXV+PMmTOoqamZ1ebrr79G\nY2MjGhoa8N577+GZZ55Zls8mIqLFWTA0Dh8+jOHhYYyNjSE+Ph4xMTH47W9/u+QPLikpQXh4OIKD\ng2FnZ4dDhw7hiy++mNXmyy+/xLFjxwAAKSkpGBwcRHd395I/m4iIFmfBA+HV1dVwd3fHH/7wB6Sl\npeHVV19FYmIi/u7v/m5JH9ze3g6ZTGZ4LJVKUVxcvGCbtrY2+Pn5zWr3+eefG66HlZSUhKSkpCXV\nRkS02pSXl6O8vHzJ77NgaMzMzGB6ehqff/45/vqv/xp2dnbLcr9wY99DEIQFX5eRkYG4uLgl10RE\ntFrd+YX6vffeW9T7LDg8dfz4cQQHB2N0dBRbt25Fc3MzPDw8FvVht5NIJGhtbTU8bm1thVQq/cE2\nbW1tvFgiEZEZLRgazz33HNrb2/HNN9/AxsYGQUFByMrKWvIHJycno6GhAc3NzZiamsInn3yC9PT0\nWW3S09Nx+vRpAEBRURE8PT3vGpoiIiLTWXB46uWXXzbcwQ/4/+Ghl156aWkfLBbj7bffRmpqKnQ6\nHZ566inExMTg3XffBfDdHs7u3bvx9ddfIzw8HC4uLnj//feX9JlERLQ0C4aGi4uLISi0Wi3Onj0L\nuVy+LB+elpaGtLS0Wc8dP3581uO33357WT6LiIiWbsHQuPOihM8//zwefPDBFSuIiIgs1z3PCB8b\nG0N7e/tK1EJERBZuwT2N+Ph4w5/1ej16enqWfDyDiIis04Kh8dVXX/1/Y7EYfn5+sLOzW9GiiIjI\nMs0bGsPDw3B3d4e7u/us50dGRgAA3t7eK1sZERFZnHlD4/Dhwzh37hwSExPnnIXd1NS0ooUREZHl\nmTc0zp07BwBobm42VS1ERGTh5g2NioqKH3xhYmLishdDRESWbd7Q+OUvfwmRSAStVovy8nIoFAoA\nQGVlJZKTk6FWq01WJBERWYZ552lkZ2cjKysLgYGBqKioMFxWV6PRIDAw0JQ1EhGRhVhwcl9tbe2s\nuRpxcXF33WGPiIjWhgXnaSgUCjz99NM4evQoBEHARx99BKVSaYraiIjIwiwYGu+//z7eeecdvPnm\nmwCArVu38l7dRERrlEi489Z4VkYkEuHUqVO8cx8R0T1ITk6+686oxlhwT6O+vh7/+I//iOrqami1\nWgDf/aK+cePGvVdJRERWbcED4U888QT+8i//EmKxGFlZWTh27BiOHDliitqIiMjCLBgaWq0WDzzw\nAARBQHBwMP75n//ZMFuciIjWlgWHpxwdHaHT6RAeHo63334bgYGBGBsbM0VtRERkYRYMjTfeeAPj\n4+M4efIk/umf/gnDw8P44IMPTFEbERFZmAWHpzZt2gQ3NzfIZDKcOnUKf/rTn9DS0mKK2oiIyMLM\nGxqjo6P43e9+h7/6q7/C73//e+j1enz22WeIjY3FH/7wB1PWSEREFmLe4ak///M/h7u7OzZv3owL\nFy7g1KlTcHR0xEcffQSVSmXKGomIyELMGxqNjY2orKwEADz99NMICAhAS0sLnJycTFYcERFZlnmH\np2xtbWf9WSKRMDCIiNa4efc0Kisr4ebmZnis1WoNj0UiEYaHh1e+OiIisijzhoZOpzNlHUREZAUW\nPOWWiIjoewwNIiIyGkODiIiMxtAgIiKjMTSIiMhoDA0iIjIaQ4OIiIzG0CAiIqMxNIiIyGgMDSIi\nMhpDg4iIjMbQICIiozE0iIjIaAwNIiIyGkODiIiMxtAgIiKjMTSIiMhoDA0iIjIaQ4OIiIzG0CAi\nIqMxNIiIyGhic3xof38/Dh48iJaWFgQHB+O///u/4enpeVe74OBguLu7w9bWFnZ2digpKTFDtURE\n9D2z7Gm8+uqr2LVrF+rr67Fz5068+uqrc7YTiUTIzs6GRqNhYBARWQCzhMaXX36JY8eOAQCOHTuG\nzz//fN62giCYqiwiIlqAWYanuru74efnBwDw8/NDd3f3nO1EIhEeeOAB2Nra4vjx4/jZz342Z7vP\nP/8chYWFAICkpCQkJSWtTOFERFaqvLwc5eXlS36fFQuNXbt2oaur667nf/Ob38x6LBKJIBKJ5nyP\ngoICBAQE4NatW9i1axeio6OxZcuWu9plZGQgLi5ueQonIlqF7vxC/d577y3qfVYsNDIzM+dd5+fn\nh66uLvj7+6OzsxPr16+fs11AQAAAwNfXF3v37kVJScmcoUFERKZhlmMa6enp+OCDDwAAH3zwATIy\nMu5qMz4+jpGREQDA2NgYLly4gPj4eJPWSUREs5klNF544QVkZmYiMjISly5dwgsvvAAA6OjowJ49\newAAXV1d2LJlC1QqFVJSUvDwww/jwQcfNEe5RET0f0SClZ+eJBKJcOrUKR7TICK6B8nJyYs6O5Uz\nwomIyGhmOeWW1hadToehoSEMDg5iZGQE4+PjmJqawszMDPR6veEMOrFYDLFYDDs7O9jb28PR0RGO\njo5wdnaGi4sLnJycYGPD7zlE5sTQoGUxMTGBpqYm3LhxAy0tLejs7ERHRwe6u7vR29sLd3d3+Pj4\nwNPTE25ubnB0dISdnZ0hBHQ6HWZmZjA9PY2pqSlotVpotVqMj49jdHQUo6OjmJiYgKurK9zd3eHp\n6Tlr8fLywrp167Bu3TqsX78e69evh7Ozs5n/VohWH4YG3bPJyUnU1NTg2rVrqK2tRV1dHTo6OhAW\nFga5XA65XI4dO3YgNDQUMpkM/v7+sLOzW/LnzszMYHh4GAMDA+jr68OtW7fQ29uL7u5udHV1oaqq\nyhBWHR0dsLOzQ2BgIPz9/eHv74+AgABIpVJIJBJIpVI4Ojouw98G0drC0KAFjY2N4fLly6ioqMCV\nK1dQV1eHqKgo3HfffThw4ACSkpIQHR29LMHwQ8RiMby9veHt7Y2wsLAfbCsIAvr7+9HS0mJYGhsb\nkZmZievXr+PmzZvw9vZGcHAwZDIZgoKCEBwcjJCQEPj6+s474ZRorWNo0F10Oh1qamqgVqtRWlqK\n2tpaJCUl4Sc/+QmOHTuGTZs2wcXFxdxl/iCRSAQfHx/4+PggMTHxrvU6nQ4tLS2or69HXV0dampq\n8PHHH6O2thaTk5MICwtDaGgoQkNDER4ejvDwcHh4eJihJ0SWhafcEgBgZGQEhYWFhiUgIAAPPfQQ\nHnroIfzoRz9aU8cH+vr6UFVVhWvXruHy5cuorKzEtWvX4ObmhqioKISFhSE6OhpRUVEICAjgXglZ\npcWecsvQWMO6urqQk5ODvLw8VFVVYcuWLUhPT8eePXsgk8nMXZ5F0ev1aG5uxuXLl6HRaFBWVgaN\nRoPJyUnExMQgKioKMTExkMvl8Pf3Z5CQxWNoMDSM0tbWhm+//RY5OTloa2vDnj17sG/fPuzatcvi\nh5wsUVdXF8rLy1FWVobi4mKUlZVBEATExsYiJiYGsbGxiI2Nhbu7u7lLJZqFocHQmNf3QXHp0iX0\n9PQgIyMDjz32GLZv377iB6/XGkEQ0NbWhpKSEhQXF0OtVkOj0cDPzw9xcXGIjY2FQqFAaGgobG1t\nzV0urWGLDQ0eCF+luru7kZmZiYsXL6Krqwv79u3DW2+9ha1bt/KX1QoSiUSQyWSQyWR49NFHAXx3\nqvC1a9dQVFSEgoICfPrpp+ju7kZcXBzi4uKgUCgQHx8PV1dXM1dPtDDuaawig4ODuHjxIi5evIjr\n168jIyPamXcMAAAPGklEQVQDhw8fxo4dOyAW8/uBJent7TWESH5+PioqKrBhwwYoFAoolUoolUr4\n+/ubu0xaxTg8tUZDY2xsDDk5OcjMzMTly5eRlpaGI0eOIDU1FQ4ODuYuj4w0NTUFjUaD/Px85OTk\nQK1Ww97eHiqVCkqlEgkJCQgJCeEBdlo2DI01FBpTU1NQq9XIzMxEQUEBfvzjH+Po0aNIT0/nEMcq\nIQgC6uvrDSGSl5eHoaEhJCQkGEIkMjKSe5C0aAyNVR4aOp0OFRUVyMzMxKVLlxAbG4sjR47gwIED\n8PHxMXd5ZALt7e3Izc1FTk4OsrOz0dHRAZVKBYVCgYSEBMTGxvLEBjIaQ2MVhoYgCLh27RouXLiA\nixcvIiAgAEeOHMGhQ4c4j4LQ29uLvLw85OTkICsrC42NjYiPj4dSqURSUhLi4uI4REnzYmisktAQ\nBAENDQ2GM58cHBzw+OOP4/Dhw4iOjjZ3eWTBhoaGkJ+fj+zsbGRnZ6O6uhpyuRwqlQoJCQlQKBRw\ncnIyd5lkIXjKrRUTBAHXr1/HxYsXcenSJczMzODgwYP44osvoFKpePCTjOLh4YE9e/YYbpk8OjqK\nwsJCZGVl4fTp07h69SoiIyOhUqmQmJgIpVLJCZ10z7inYSbfH+jMyspCVlYWJicn8dhjj+HgwYPY\ntGkTg4KW3fj4ONRqNbKzs5GVlQWNRoPw8HBDiKhUKp5IsYZweMoKQkOn0+Hq1auGA5k2NjZ47LHH\nsH//fgYFmZxWq0VxcbEhRMrKyhASEjIrRHhl39WLoWGhoTExMYHi4mLk5+cjLy8P/v7+yMjIwP79\n+6FQKBgUZDEmJydRWlpqOLBeXFwMiURiOCaiUqmwbt06c5dJy4ShYUGh0dbWhoKCAhQVFUGj0SAx\nMRF79+5Feno6QkNDzV0ekVGmp6eh0WgMB9YLCwvh7e1tmCuSmJiIgIAAc5dJi8TQMGNojI6Oory8\nHKWlpSgqKsL4+DhSU1PxyCOP4MEHH+QuPq0KOp0OVVVVyMnJMUw4tLOzM8xaVyqVCA0NNdz3nSwb\nQ8OEoTE+Po4rV66goqICFRUVaGxsxMaNG/HQQw8hNTUVCoWC/3Fo1fv+9PDv54rk5eVhcHAQKpUK\n8fHxUKlUiImJgb29vblLpTkwNFYoNARBQFdXFyorK1FVVYWrV6/i+vXrSExMxI4dO7Bjxw7cf//9\ncHR0XJHPJ7ImHR0dKCgoMIRIY2MjYmJiEB8fD4VCAYVCAU9PT3OXSWBoLEto6PV6dHR0oK6uDvX1\n9aivr0d1dTVsbW2RkpKC+++/H1u2bEFSUhJDgsgIIyMjKC4uRl5eHvLz81FaWor169cjPj7esAQH\nB3PP3AwYGkaGhiAIGBoaQnt7O9ra2tDW1oabN2+ipaUFN27cgLe3NxQKBRITE5GUlISNGzciMDCQ\nZzkRLYPvTzsvKChAQUEB1Go1+vv7oVAoIJfLER8fzzsdmsiaDo3nn38eMpkM09PTmJychFarxdjY\nGEZHRzE8PIzBwUH09fXh1q1b6OrqgqOjI2QyGcLDwxEZGWm4x7NcLuc/ViIT6+7uhlqtRmFhoeFO\nh/7+/oiNjYVcLkdcXBzCw8N5Rd9ltqZDY/v27bC1tYW9vT2cnJzg4uICDw8PeHl5wdfXFz4+PvD3\n90dAQACkUinc3NzMXTYRzWNmZgZVVVUoLi5GYWEhSkpK0NLSgujoaMTExBiWDRs2cFhrCdZ0aFh5\nF4hoAcPDw6ioqEBxcTFKSkpQVlaGgYEByOVyREZGGgJFKpUySIzE0CCiNaW3t9cwP6q0tBQVFRUY\nGhpCdHQ0IiIiEBkZiaioKISEhHBoaw4MDSJa8/r6+gzzp8rLy3HlyhXcvHkTISEhiIiIQFhYGMLD\nwxEREQFvb29zl2tWDA0iojmMj4+jqqoKV65cgUajMcy5srOzQ0REBEJCQhAaGmpY1sqVfhkaRERG\nEgQB7e3tuHr1qmHSblVVFerq6uDu7o6QkBAEBQUhKCgIwcHBCAoKgq+v76o69Z6hQUS0RHq9Hjdv\n3kRNTQ2qq6tRW1uL6upqNDQ0YHx8HEFBQZDJZJBKpYafUqkUPj4+VhcoDA0iohU0ODiI+vp6NDQ0\noKGhwfDnpqYmjI+PQyqVQiKRGE7vv/2nl5eXxYUKQ4OIyEyGh4fR1NRkWG7cuIGmpibcvHkTra2t\n0Gq18Pf3h5+fH3x9fbFu3TqsW7fOMI/Mx8cH3t7ecHFxMVm4MDSIiCzU2NgY2tra0N7ejtbWVrS3\ntxuWrq4udHd3o6enBzqdDj4+PvDy8oKXlxc8PDzg5uYGd3d3uLu7w83NDa6urnB1dYWLiwtcXV3h\n5OQEZ2dnODg43NMcFYYGEZGVGxsbw61bt3Dr1i309vait7cXfX196O/vR29vL4aGhjAwMIChoSEM\nDQ1hdHQUY2NjGBsbw8TEBBwcHODo6AgHBwc4ODjA3t4eYrEYYrEYtra2sLW1hY2NDWxsbFBSUsLQ\nICJaq/R6PSYmJqDVaqHVajExMYHJyUlMT09jenoaMzMz0Ol00Ol0EAQBO3bsYGgQEZFxFvu7kxdp\nISIiozE0iIjIaAwNIiIyGkODiIiMxtAgIiKjMTQsXHZ2trlLWDGruW8A+2ftVnv/FsssofHpp58i\nNjYWtra2qKiomLfd+fPnDTdUee2110xYoeVYzf9wV3PfAPbP2q32/i2WWUIjPj4en332GbZu3Tpv\nG51Oh2effRbnz59HdXU1zpw5g5qaGhNWSUREdzLLPRCjo6MXbFNSUoLw8HAEBwcDAA4dOoQvvvgC\nMTExK1wdERHNSzCj7du3C+Xl5XOu+/TTT4Wnn37a8PjDDz8Unn322bvaAeDChQsXLotYFmPF9jR2\n7dqFrq6uu55/5ZVX8Mgjjyz4emMvDyzwEiJERCazYqGRmZm5pNdLJBK0trYaHre2tkIqlS61LCIi\nWgKzn3I7355CcnIyGhoa0NzcjKmpKXzyySdIT083cXVERHQ7s4TGZ599BplMhqKiIuzZswdpaWkA\ngI6ODuzZswcAIBaL8fbbbyM1NRVyuRwHDx7kQXAiIjMzS2js3bvXcAvErq4ufPPNNwCAwMBAnDt3\nztAuLS0NdXV1aGxsREJCwoJzNp577jlERERAqVRCo9GYpC/LZaE5KdnZ2fDw8EBCQgISEhJw4sQJ\nM1S5OE8++ST8/PwQHx8/bxtr3nYL9c+at11rayt27NiB2NhYxMXF4eTJk3O2s9btZ0z/rHn7TUxM\nICUlBSqVCnK5HP/wD/8wZ7t72n6LOnxuYjMzM0JYWJjQ1NQkTE1NCUqlUqiurp7V5ty5c0JaWpog\nCIJQVFQkpKSkmKPURTGmf1lZWcIjjzxipgqXJjc3V6ioqBDi4uLmXG/N204QFu6fNW+7zs5OQaPR\nCIIgCCMjI0JkZOSq+r9nTP+sefsJgiCMjY0JgiAI09PTQkpKipCXlzdr/b1uP7Mf0zDG7XM27Ozs\nDHM2bvfll1/i2LFjAICUlBQMDg6iu7vbHOXeM2P6B1jvmWJbtmyBl5fXvOutedsBC/cPsN5t5+/v\nD5VKBQBwdXVFTEwMOjo6ZrWx5u1nTP8A691+AODs7AwAmJqagk6ng7e396z197r9rCI02tvbIZPJ\nDI+lUina29sXbNPW1mayGpfCmP6JRCIUFhZCqVRi9+7dqK6uNnWZK8aat50xVsu2a25uhkajQUpK\nyqznV8v2m69/1r799Ho9VCoV/Pz8sGPHDsjl8lnr73X7mWVG+L1a7JwNY19nbsbUmZiYiNbWVjg7\nO+Obb75BRkYG6uvrTVCdaVjrtjPGath2o6Oj2L9/P9588024urretd7at98P9c/at5+NjQ0uX76M\noaEhpKamIjs7G9u3b5/V5l62n1XsaRgzZ+PONm1tbZBIJCarcSmM6Z+bm5thNzMtLQ3T09Po7+83\naZ0rxZq3nTGsfdtNT0/j0UcfxdGjR5GRkXHXemvffgv1z9q33/c8PDywZ88elJWVzXr+XrefVYSG\nMXM20tPTcfr0aQBAUVERPD094efnZ45y75kx/evu7jZ8GygpKYEgCHeNTVora952xrDmbScIAp56\n6inI5XL84he/mLONNW8/Y/pnzduvt7cXg4ODAACtVovMzEwkJCTManOv288qhqdun7Oh0+nw1FNP\nISYmBu+++y4A4Pjx49i9eze+/vprhIeHw8XFBe+//76ZqzaeMf374x//iHfeeQdisRjOzs74+OOP\nzVy18Q4fPoycnBz09vZCJpPh5ZdfxvT0NADr33bAwv2z5m1XUFCA//qv/4JCoTD8snnllVdw8+ZN\nANa//YzpnzVvv87OThw7dgx6vR56vR4//elPsXPnziX97hQJ1nxaABERmZRVDE8REZFlYGgQEZHR\nGBpERGQ0hgYRERmNoUG0zCYnJ7Ft27Z7uvTEyZMn8eGHH65gVUTLg2dPES2z//zP/0RfXx+ef/55\no18zMjKCnTt3oqSkZAUrI1o67mkQGam0tBRKpRKTk5MYGxtDXFzcnNchOnPmDP7sz/4MwHeX1d62\nbRsyMjIQFhaGF154AR9++CE2bdoEhUKBGzduAPhu1rGPjw+uXbtm0j4R3SurmNxHZAk2btyI9PR0\nvPjii9BqtfjpT39618XfdDodqqqqEBkZaXiusrIStbW18PLyQkhICH72s5+hpKQEJ0+exFtvvYV/\n//d/BwBs2rQJubm5iI2NNWm/iO4FQ4PoHrz00ktITk6Gk5MT3nrrrbvW9/b2ws3NbdZzGzduNFyW\nITw8HKmpqQCAuLg4ZGVlGdoFBgYa9jyILBWHp4juQW9vL8bGxjA6OgqtVjtnmzsPEzo4OBj+bGNj\nY3hsY2ODmZmZWa+ztqvD0trD0CC6B8ePH8eJEyfw+OOP4+///u/vWr9u3TqMjo4u6r07OzsRHBy8\nxAqJVhZDg8hIp0+fhoODAw4dOoQXXngBpaWlyM7OntXG1tYWcXFxqKurA/DdfQnm23u4c11JSQm2\nbNmyYvUTLQeecku0zE6dOoXu7u4590TmMzw8jJ07d6K0tHQFKyNaOoYG0TKbmprCAw88gJycHKOP\nUZw8eRLe3t44evToCldHtDQMDSIiMhqPaRARkdEYGkREZDSGBhERGY2hQURERmNoEBGR0RgaRERk\ntP8FlAv6KSou0D8AAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 1. Back pressure pB is ambient pressure\n",
      "p_B = p_amb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 2. Find minimum area in domain\n",
      "A_e = A(x[-1])  # m^2\n",
      "A_c = np.min(A(x))  # m^2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 3. Compute subsonic and supersonic Mach number at the exit\n",
      "fl = isentropic.IsentropicFlow(gamma=g)\n",
      "M_e_sub, M_e_sup = isentropic.mach_from_area_ratio(fl, A_e / A_c)\n",
      "\n",
      "print(A_e / A_c)\n",
      "print(M_e_sub, M_e_sup)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2.5\n",
        "0.23954284305818874 2.44276484552\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/home/juanlu/.local/lib/python3.3/site-packages/skaero/gasdynamics/isentropic.py:150: RuntimeWarning: divide by zero encountered in double_scalars\n",
        "  ((self.gamma + 1) / (2 * (self.gamma - 1))) / M\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 4. Compute pressure limit values\n",
      "from IPython.core.display import Latex\n",
      "\n",
      "Latex(r\"\\[\\frac{p_0}{p} = (1 + \\frac{\\gamma - 1}{2} M^2)^{\\gamma / (\\gamma - 1)}\\]\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "\\[\\frac{p_0}{p} = (1 + \\frac{\\gamma - 1}{2} M^2)^{\\gamma / (\\gamma - 1)}\\]"
       ],
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<IPython.core.display.Latex at 0x7fe7d808b390>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Isentropic limit subsonic and supersonic solutions\n",
      "p_e3 = p_0 * fl.p_p0(M_e_sub)\n",
      "p_e6 = p_0  * fl.p_p0(M_e_sup)\n",
      "\n",
      "print(\"Exit pressure of the limit isentropic solutions\")\n",
      "print(p_e3, p_e6)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Exit pressure of the limit isentropic solutions\n",
        "194716.088639 12966.4263803\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Normal shock at the exit of the nozzle\n",
      "# ..."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Isentropic pressure solutions\n",
      "M_isen = np.empty((len(x), 2))\n",
      "for i in range(len(x)):\n",
      "    M_isen[i] = isentropic.mach_from_area_ratio(fl, A(x[i]) / A_c)\n",
      "\n",
      "# Subsonic isentropic limit solution\n",
      "M_isen_sub = M_isen[np.arange(len(x)), np.zeros(len(x), dtype=int)]\n",
      "p_isen_sub = fl.p_p0(M_isen_sub)\n",
      "\n",
      "# Supersonic isentropic solution\n",
      "i_c = np.argmin(A(x))\n",
      "mask = np.zeros(len(x), dtype=int)\n",
      "mask[i_c:] = 1\n",
      "M_isen_sup = M_isen[np.arange(len(x)), mask]\n",
      "p_isen_sup = fl.p_p0(M_isen_sup)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 5. Discriminate flow\n",
      "# TODO: I need to compute also the shock inside nozzle, overexpanded, properly expanded\n",
      "# and underexpanded cases.\n",
      "if p_B > p_e3:\n",
      "    print(\"Fully subsonic flow\")\n",
      "elif p_B <= p_e3:\n",
      "    print(\"Chocked flow\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Chocked flow\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Given normal shock location, return pressure\n",
      "def pe_p0_from_shock_location(A_loc_ratio, A_e_ratio):\n",
      "    # 1. Compute M_1\n",
      "    M_1 = isentropic.mach_from_area_ratio(fl, A_loc_ratio)[1]\n",
      "    # 2. Compute M_2\n",
      "    ns = shocks.NormalShock(M_1, gamma)\n",
      "    M_2 = ns.M_2\n",
      "    # 3. Compute A_2 over A_2^*\n",
      "    a2_a2c = fl.A_Astar(M_2)\n",
      "    # 4. Compute A_e / A_2^*\n",
      "    ae_a2c = A_e_ratio / A_loc_ratio * a2_a2c\n",
      "    # 5. Compute M_e\n",
      "    M_e = isentropic.mach_from_area_ratio(fl, ae_a2c)[0]\n",
      "    # 6. Compute p_e over p_0e\n",
      "    p0e_pe = 1 / fl.p_p0(M_e)\n",
      "    # 7. Compute p_02 over p_2\n",
      "    p02_p2 = 1 / fl.p_p0(M_2)\n",
      "    # 8. Compute p_2 over p_1\n",
      "    p2_p1 = ns.p2_p1\n",
      "    # 9. Compute p_01 over p_1\n",
      "    p01_p1 = 1 / fl.p_p0(M_1)\n",
      "    # 10. Compute p_02 over p_01\n",
      "    p02_p01 = p02_p2 * p2_p1 * 1 / p01_p1\n",
      "    # 11. Compute p_e over p_0\n",
      "    pe_p0 = 1 / p0e_pe * p02_p01\n",
      "    return pe_p0\n",
      "\n",
      "\n",
      "#print(p_e_from_shock_location(2.0, 3.0))  # Anderson example 5.6, it works\n",
      "def F(A_loc_ratio, A_e_ratio, pe_p0):\n",
      "    return pe_p0_from_shock_location(A_loc_ratio, A_e_ratio) - pe_p0\n",
      "\n",
      "#print(sp.optimize.brentq(F, 1.0, 3.0, args=(3.0, 0.5)))  # Anderson example 5.7, works\n",
      "#print(F(2.3397, 3.0, 0.5))  # 5.84275393312e-06\n",
      "\n",
      "A_shock = sp.optimize.brentq(F, 1.0, 3.0, args=(A_e / A_c, p_B / p_0))\n",
      "print(\"Location of the normal shock\")\n",
      "print(A_shock * A_c)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Location of the normal shock\n",
        "0.0222720756212\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Find corresponding index\n",
      "def find_nearest_index(a, a0):\n",
      "    \"\"\"Element in nd array a closest to the scalar value a0.\n",
      "\n",
      "    See http://stackoverflow.com/a/10465997/554319\n",
      "\n",
      "    \"\"\"\n",
      "    idx = np.abs(a - a0).argmin()\n",
      "    return idx\n",
      "\n",
      "\n",
      "Area = A(x)\n",
      "i_s = find_nearest_index(Area, A_shock * A_c)\n",
      "x[i_s]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "2.46"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute pressure along nozzle axis\n",
      "p_ratio = np.empty_like(x)\n",
      "\n",
      "# Isentropic portion\n",
      "p_ratio[:i_s] = p_isen_sup[:i_s]\n",
      "\n",
      "# Portion behind the shock\n",
      "# 1. Compute M_1\n",
      "M_1 = isentropic.mach_from_area_ratio(fl, A_shock)[1]\n",
      "# 2. Compute M_2\n",
      "ns = shocks.NormalShock(M_1, gamma)\n",
      "M_2 = ns.M_2\n",
      "# 3. Compute A_2 over A_2^*\n",
      "a2_a2c = fl.A_Astar(M_2)\n",
      "A_c2 = A_shock * A_c / a2_a2c\n",
      "\n",
      "M_post = np.empty_like(x[i_s:])\n",
      "for i in range(len(x[i_s:])):\n",
      "    M_post[i] = isentropic.mach_from_area_ratio(fl, A(x[i_s + i]) / A_c2)[0]\n",
      "\n",
      "print(M_post)\n",
      "p_ratio[i_s:] = fl.p_p0(M_post)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.53107106  0.52502108  0.51928679  0.51385348  0.50870794  0.50383834\n",
        "  0.49923405  0.49488553  0.49078421  0.48692242  0.48329332  0.47989083\n",
        "  0.47670958  0.47374485  0.47099257  0.46844926  0.46611202  0.46397849\n",
        "  0.46204687  0.4603159   0.45878484  0.45745348  0.45632216  0.45539176\n",
        "  0.45466371  0.45414004  0.45382338  0.45371699]\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Mach number distribution\n",
      "M = np.zeros_like(x)\n",
      "M[:i_s] = M_isen_sup[:i_s]\n",
      "M[i_s:] = M_post"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot whole solution\n",
      "\n",
      "fig = plt.figure(figsize=(6, 8))\n",
      "\n",
      "#ax_nozzle = fig.add_subplot(211)\n",
      "#ax_nozzle.fill_between(x, radius, -radius, facecolor=\"#cccccc\")\n",
      "\n",
      "ax_press = fig.add_subplot(211)\n",
      "ax_press.plot(x[i_c:], p_isen_sub[i_c:], 'k--', x[i_s:], p_isen_sup[i_s:], 'k--')\n",
      "ax_press.annotate(s=\"Subsonic isentropic\\nsolution\", xy=(2.0, 0.9), xytext=(1.5, 0.7), arrowprops=dict(arrowstyle = \"->\"))\n",
      "ax_press.annotate(s=\"Supersonic isentropic\\nsolution\", xy=(1.5, 0.2), xytext=(0.1, 0.1), arrowprops=dict(arrowstyle = \"->\"))\n",
      "ax_press.set_ylabel(r\"$p / p_0$\")\n",
      "ax_press.plot(x, p_ratio)\n",
      "ax_press.set_title(\"Pressure\")\n",
      "\n",
      "ax_mach = fig.add_subplot(212)\n",
      "ax_mach.plot(x, M)\n",
      "ax_mach.set_ylabel(\"$M$\")\n",
      "ax_mach.set_title(\"Mach number\")\n",
      "\n",
      "fig.suptitle(\"Exercise 18\")\n",
      "fig.savefig(\"exercise_18.png\", dpi=100)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAIGCAYAAABH+rhpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8zdf/wPHXzZCQRIYRWcSIhMpErCI2tVqlaGvVbq3q\nV1UH0WGVryqqqFpVpUbVqFKbIkas2iVGCBkyJDLuvef3h2/vTyQ0lyQ34/18PPJIbu65n/P+HHHf\n93zO55yjUUophBBCiH9hZuoAhBBCFA6SMIQQQuSIJAwhhBA5IglDCCFEjkjCEEIIkSOSMIQQQuSI\nJAxR5JibmxMYGGj4mjZtWp7VNWHCBHbs2PFcx4iLi6NZs2bY2dkxfPjwTM8tXrwYX19f/P39adeu\nHbGxsc9VlxDPQyPzMERRY2dnR1JSUq4d75//IhqNJteO+aiUlBTCw8M5c+YMZ86cYfbs2QCkp6fj\n4uLCpUuXcHJyYuzYsZQqVYoJEybkSRxC/BvpYYhiISEhAR8fHy5evAhAz549WbRoEQBffvklwcHB\n+Pv7ExoaCkBERATe3t706dMHX19fbty4wdSpU/Hz8yMgIIAPP/wQgL59+7J27VoAPvjgA1544QX8\n/f0ZM2YMANHR0XTt2pXg4GCCg4P5888/s8RWqlQpGjVqhJWVVabfW1hY4OjoyP3791FKkZiYiJub\nW560jxA5YWHqAITIbQ8ePCAwMNDw+MMPP6Rbt27MmTOHvn37MmLECBISEujfvz/btm3j8uXLhIWF\nodfr6dy5M/v27cPDw4PLly+zfPlygoOD+e233/j1118JCwvD2tqa+Ph44GGvQ6PREBsbyy+//ML5\n8+cBSExMBGDkyJG8++67NGrUiOvXr9O2bVvOnj2bbdyP92DMzMyYNWsWtWrVwtbWlurVqzN37ty8\naDIhckQShihySpYsSXh4eJbft2zZktWrVzNs2DBOnToFwLZt29i2bZshwSQnJ3P58mU8PDyoVKkS\nwcHBAOzYsYO33noLa2trABwcHDId28HBAWtra/r370+HDh3o0KEDAH/88Qfnzp0zlEtKSiIlJYVS\npUr963kkJiYyYsQITp48SeXKlRk+fDiTJ0/mo48+eoZWEeL5ScIQxYZer+fcuXPY2NgQFxeHq6sr\nAOPGjWPQoEGZykZERGBjY5Ppd08a7lNKYW5uTlhYGDt27GDNmjXMmTOHHTt2oJTi8OHDlChRwuh4\nz507R+XKlalcuTIA3bp1Y+rUqUYfR4jcImMYotiYOXMmL7zwAitWrKBfv35otVratGnD999/T3Jy\nMgCRkZFER0dneW2rVq1YvHgxDx48AODevXuZnk9OTiY+Pp527drx3//+l5MnTwLQunVrvv76a0O5\nEydOPDG+xxNSlSpVOH/+PDExMQBs376dmjVrPsOZC5E7pIchipzHxzDatWtH3759WbRoEUeOHMHG\nxoYmTZrwxRdfMGHCBM6dO0eDBg2Ah3dY/fDDD4axiX+0adOGEydOUKdOHUqUKEH79u35/PPPgYdj\nD0lJSXTu3JnU1FSUUsycOROAr7/+mnfeeQd/f3+0Wi1Nmzblm2++yRKzp6cnSUlJpKen88svv7B9\n+3Z8fHyYNGkSzZo1w8zMDE9PT5YsWZKHLSfE08lttUIIIXJELkkJIYTIEUkYQgghckQShhBCiByR\nhCGEECJHJGEIIYTIEUkYQgghckQShhBCiByRhCGEECJHJGEIIYTIEUkYQgghckQShhBCiByRhCGE\nECJHJGEIIYTIEUkYQgghckQShhBCiByRhCGEECJHJGEIIYTIEUkYQgghckQShhBCiByRhCGEECJH\nJGEIIYTIEUkYQgghckQShhBCiByRhCGKHU9PT0qVKoWdnR0VKlSgX79+JCcnmzosIQo8SRii2NFo\nNGzatImkpCSOHz/O0aNH+fzzzzOV0Wq1Joru/+l0OlOHIEQmkjBEsebq6kq7du04c+YMZmZmfPPN\nN3h5eeHt7Q3Apk2bCAgIwNHRkUaNGnH69GnDa6dOnYq7uzulS5fGx8eHnTt3AhAWFkadOnWwt7en\nQoUKvPfeewDs3r0bDw+PTPV7enoaXhcaGkrXrl3p1asX9vb2LF26lISEBPr374+rqyvu7u588skn\n6PX6/GgaIbKQhCGKJaUUADdu3GDLli0EBgYCsGHDBo4cOcLZs2cJDw+nf//+LFy4kLi4OAYPHkyn\nTp3IyMjgwoULzJ07l6NHj5KYmMi2bdvw9PQEYOTIkbz77rskJCRw5coVunfv/sQ4NBpNpse//vor\n3bp1IyEhgddff52+fftSokQJ/v77b8LDw9m2bRvfffdd3jSKEP9CEoYodpRSvPzyyzg6OtK4cWNC\nQkL48MMPARg3bhwODg5YWVmxYMECBg8eTN26ddFoNPTu3RsrKysOHjyIhYUFaWlp/PXXX2RkZFCx\nYkWqVKkCQIkSJbh06RIxMTGUKlWK4ODgHMfWsGFDOnXqBEBCQgK//fYbM2fOpGTJkpQrV45Ro0bx\n008/5X6jCJEDkjBEsaPRaNiwYQP37t0jIiKCOXPmYG1tDZDpktG1a9eYMWMGjo6Ohq+bN29y+/Zt\nqlatyldffUVoaCjOzs707NmT27dvA7Bo0SIuXrxIjRo1CA4OZvPmzTmOzd3dPVP9GRkZuLi4GOof\nMmQI0dHRudQSQhhHEoYQj3j0ElHFihX56KOPuHfvnuHr/v37hktMPXv2ZN++fVy7dg2NRsPYsWMB\nqFatGj/++CPR0dGMHTuWrl278uDBA2xsbEhJSTEcX6fTZXnzf7R+Dw8PrKysiI2NNdSfkJCQaRxF\niPwkCUOIJxg4cCDffvstYWFhKKVITk5m8+bN3L9/n4sXL7Jz507S0tKwsrLC2toac3NzAH744QdD\nIrC3t0ej0WBmZkb16tVJTU1ly5YtZGRk8Pnnn5OWlvbE+l1cXGjdujWjR48mKSkJvV7P33//zd69\ne/Pl/IV4nCQMIf7n8QHo2rVrs3DhQoYNG4aTkxNeXl4sW7YMgLS0NMaNG0e5cuVwcXEhJiaGyZMn\nA/D7779Tq1Yt7OzsePfdd/npp5+wsrLC3t6eb775hgEDBuDu7o6trW2mS2AajSZLDMuWLSM9PZ2a\nNWvi5OREt27diIqKyuOWECJ7GvXP7SJCCCHEU+R7D+Ott97C2dkZX1/fJ5YZMWIEXl5e+Pv7Ex4e\nno/RCSGEeJJ8Txj9+vVj69atT3x+y5YtXL58mUuXLrFgwQKGDh2aj9EJIYR4Eov8rrBx48ZEREQ8\n8flff/2VPn36AFCvXj3i4+O5c+cOzs7Omco9fq1XCCFEzjzrSESBG/SOjIzMNBDo7u7OzZs3sy2r\nlCqyXxMmTDB5DHJ+cm5yfkXv63kUuIQBWbOf9CaEEML0ClzCcHNz48aNG4bHN2/exM3NzYQRCSGE\ngAKYMDp16mS41/3QoUM4ODhkGb8oDkJCQkwdQp4qyudXlM8N5PyKs3yfh9GzZ0/27NlDTEwMzs7O\nTJw4kYyMDAAGDx4MwLBhw9i6dSs2NjYsXryYoKCgrIFrNM99PU4IIYqb53nvLLQT9yRhCCGE8Z7n\nvbPAXZISQghRMOX7PAwhhBC5JyUlheTkZLRaLTqdDr1ej16vx8HBAQcHh1ytSxKGEEI8h7S0NBIT\nE7GysqJ06dJZnj906BB//vknKSkphq8HDx7QqVMn2rdvn6X8V199xbx580hPTyc9PZ20tDQyMjIY\nP368YbvfR02fPp2vv/4aCwsLzM3NMTc3x8zMjLFjx+b6ShkyhiGEKNa0Wi1xcXHExcURGxtr+Llm\nzZrUrVs3S/k5c+bw3//+l8TERJKSkgCws7Pjk08+YeTIkVnKb9y4kZ07d2JjY0PJkiUpVaoUJUuW\npEGDBvj7+2cpHxUVRXx8PFZWVlhZWWFpaUmJEiUoWbIkJUqUeO7zlUFvIYT4H6UU8fHxREVFERUV\nxZ07d7hz5w61atWiRYsWWcpPmTKF6dOnU6ZMGcqUKYOTkxNOTk507drVsF3uo+7cucP9+/ext7fH\nzs4OKyur/DitXCMJQwhRLOh0Om7fvs3169exs7PLdtXrr7/+mk8++QQXFxecnZ2pUKEC5cuX56WX\nXqJdu3YmiLpgKbYJY+Oes1hYmGFpYY6luRkWFuZY/u+xXSkrSttaUdrGGqsSMlQjRGG1c+dOJk6c\nyLVr17h9+zZlypTBw8ODN954gxEjRmQpr9frMTOTG0CfpNgmjLcnrUer06HV6knXPvyeodWRodWR\nlJJOUnIqCfdT0Wg0lLaxprStFQ52JXF2ssW5jB0VythS3smWCmXtcC5jh2tZO6ytLE19akIUeTEx\nMZw/f54LFy5w8eJFLl68SJUqVZgxY0aWsjdv3uTixYtUqlQJd3f3QncJqKAptgkjJ6ErpUhL15KY\nnEbC/VTuJT7gbtx97sTd505sElExSYafb8ckUc7RlspujlRxc6KKmxOV//flVt4eMzNZBFGI57V3\n7146deqEj48P3t7eeHt7U716dXx9ffH29jZ1eEWeJIxcotXpibybwJWbcVyNfPh1JTKOKzfjSH6Q\nTo3K5alZ1ZkXqjrjW60CXhXLShIRgoeXgS5fvkx4eDgnTpzgxIkTpKWlsXPnzmzLZrd/ucgfkjDy\nwb3EB5y9coe//r7D2St3OHXxNjEJKQR6u1Knpju1a7oT5OMql7REsZOQkEDFihVxcnIiMDCQwMBA\nAgIC8PPzo1KlSqYOTzxGEoaJxCakcPzcTY6ejeTo2ZtciIgm0NuVxkGVeTHQEx/P8tIDEYVeeno6\n4eHh7N+/n+HDh2c7FyAuLg4nJycTRCeMJQmjgEhMTuXw6evsD49g3/GrJKdm0Lq+F20aelPP1wNL\nC3NThyhEjhw6dIht27axZ88ewsLCqFq1Ki+++CKfffYZjo6Opg5PPAdJGAXUlcg4tv15kd8PXuTa\n7Xs0r1uNl5u9QAO/ipiby21/ouAaPXo0FhYWNG3alEaNGuX6mkTCdCRhFAK3YxLZeuAC63aeITY+\nhS7Na9GlpS9V3KQbL/JfQkIC27dvx93dnfr165s6HJGPJGEUMucj7rJuxxnW7/yLKu5O9Otch1b1\nvKTXIfLU9evX2bBhAxs2bCAsLIxGjRoxevRoWrVqZerQRD6ShFFIZWh1/P7nRb7fcIToe8n07VSH\nbq18KW1jberQRBHz+++/8+abb9KhQwc6d+5Mq1atsLGxMXVYwgQkYRQB4ecjWbzhKPvCr/Jm+yAG\nvByMvZ0kDpE7MjIy0Gg0WFg83zI5YWFhDBw4kJMnT+ZSZCK/ScIoQm7eSWDOTwfYdugSfTvWpm/n\nOtLjEP/q/v37rF27llWrVrFq1Srs7OxyvY7o6Gj8/Py4c+cOUVFRlC9fPtfrEHlPtmgtQtyd7Zky\n8iXWzejNtah4mg9cwMJ1YWRodaYOTRQwSikOHz7MgAED8PDw4Oeff6Zfv355stZSamoqL7/8Mu7u\n7lSrVo1Dhw7leh2i4JMeRgH3941YPlu4g8i7CUwY3IoXAz1NHZIoIMaMGcP69esZMGAAvXv3xtXV\nNc/q6tWrF2lpaZw/f5569epRtmxZJk+enGf1iZzL0OpIfpDO/ZR0klLSSEvXkpahJSNDT3qGlnSt\njoyMh4uy6hV0a+Unl6SKMqUUOw5f5tOFf1CragU+HNAc9/L2pg5LmFh8fDylS5fOl6W8p0+fTo8e\nPXjhhRdYuXIl06ZNY/fu3Xleb3GUlJLG7ehEomKSiElIIS4hhdj4/33/31dC0gOSUtK5/yANrVaP\nbakS2JaywrZkCaytLClhYY6lpRklLC2wtPjfd3MzNGYa/vteR0kYxUFqWgYL1h5mycZjjH6zMW+8\nFCgLuBVxSinOnTtHzZo1TR0K69evZ+HChaxdu5bVq1fTp08fU4dUKKWmZRBx+x5XI+8RcSuOyLuJ\n3IpO5HZ0IpHRiWi1elzK2eFS1o6yDjaUsS+Fk30pyjjYUMa+JGXsbXCws8bOxgrbklZYW1kY9T4g\ng97FzN83Yhk9YyMOdiWZOvIlKpTN/QFOYVo6nY5169YxZcoU0tPTOXbsWK7s5/w8Ro4ciZubG++/\n/75J4ygsEpJSOR9xl/NX73LpRixXI+OIuHWP2IQUKlZwwNPVEU9XRzwqOOBS1g7XcqVxLVcae1vr\nPP0gKAmjGMrQ6pi3+iDLNh1n/OCWdGpq+k+g4vllZGSwYsUKJk+ejJOTE+PGjaNDhw4FYgc5f39/\nFixYQL169UwdSoGilCIyOpGTF25z/updzl19mCTi76dSvVJZalQuT/WKZQ1767iWK23SSbqSMIqx\nU5du896MTfhXd+Xzd1rL8uqF3EcffcTBgwf55JNPCAkJKTCXHGNjY6lcuTKxsbFYWhbvv7HUtAxO\nX44i/Pwtws9Hcvz8LfR6RYC3KzWrlMencnlqVC5PxQoOBXK1akkYxVxKajofzPqNq7fimPdRFxkQ\nL8TS09NNfukpO+vXr2fBggX89ttvpg4l36VlaDl+LpL94REcOBHBxWsxVK9UlgBvV4JquBHo7Yq7\ns32BSe7/RhKGQCnFd+uPsGDtYWa934mG/rJxjcg9I0eOxNXVlbFjx5o6lDyn1yvOXb3LgRMR7D8R\nwfFzkVSvVJZGAZ408vckwNulUPfkJWEIgwMnIhj15UaGdKvPW53rFJpPPcXJuXPnGDNmDOPHjyc4\nONjU4eRIUR+/SEvXcuBEBNsPX2LH4cvYlrKicaAnjQI8qe9bkdK2RWe1BUkYIpObdxLoP/FnmgRV\nZtxbzQvkddTiKD4+ngkTJvDjjz8ybtw4hg0bViAvPz2uKI9fKKUYO2sLWw9cpEbl8rSs70Wr+l54\nuhbdTaJkaRCRibuzPaumvkH4+VuM+WqzLCtiYnq9nu+//x4fHx9SU1M5e/Yso0ePLhTJAmDv3r00\natSoyCULgHStjvU7/2L3d4NZNe0NBnYJLtLJ4nmZJGFs3boVHx8fvLy8mDp1apbnY2JiaNu2LQEB\nAdSqVYslS5bkf5CFnINdSZZ/3oO4hBSGfrGe1LQMU4dUbCUmJrJ27Vo2b97M/PnzKVeunKlDMopS\nil69epk6jDyh0+mxsDDDyb6UqUMpFPL9kpROp8Pb25s//vgDNzc36taty8qVK6lRo4ahTGhoKGlp\naUyePJmYmBi8vb25c+dOpqWZ5ZJUzmRodbz/1RYi7ybw3fiuReparBDPKzE5lUZ9vuH0mtGmDiXf\nFKpLUmFhYVSrVg1PT08sLS3p0aMHGzZsyFTGxcWFxMRE4OGnszJlyjz3Ov7FlaWFOTNGd6BGFWd6\nf7KKxORUU4ckRIGh1yvZ6dII+f4uHBkZiYeHh+Gxu7s7hw8fzlRm4MCBNG/eHFdXV5KSkli9enW2\nxwoNDTX8HBISQkhISF6EXOiZmWkIHdySCfO202/Czyz99DVsS+X+EtjF3YMHD5g3bx4jRoyQDziF\nhFanx6KIJ4zdu3fn2kKR+f5XnZPbPCdNmkRAQAC7d+/m77//plWrVpw8eTLLpjCPJgzxdBqNhtAh\nrfho7lYGfLqWpZ++hlUJeVPLLceOHePNN9/E39+fBw8e5MkGRiL36XT6In8X4eMfpidOnPjMx8r3\n1Orm5saNGzcMj2/cuIG7u3umMn/++SfdunUDoGrVqlSuXJkLFy7ka5xFkZmZhs/fbkM5RxtGTNuA\nVqc3dUiFnk6nY9KkSbRr147x48fz008/GZLFF198Qa1atfD39ycwMJCwsLCnHis0NJQZM2bkWazz\n589n+fLlOS7fqFGjXK1/z549HDx4MNeOt3HjxmxvmjFGcehh5KZ8/4hZp04dLl26REREBK6urqxa\ntYqVK1dmKuPj48Mff/xBo0aNuHPnDhcuXKBKlSr5HWqRZG5uxozRHeg/8Wc+nvs7k4e3lcl9zygx\nMZHOnTujlOLYsWOZLrUePHiQzZs3Ex4ejqWlJXFxcaSlpT31eHn97zB48GCjyh84cCBX69+1axd2\ndnY0aNAgy3M6nQ5zc3OjjtexY0c6duz4XDHJGIZx8r2lLCwsmDNnDm3atKFmzZp0796dGjVqMH/+\nfObPnw/Ahx9+yNGjR/H396dly5ZMmzYNJyen/A61yCphac63H3Xhr7+jmPezbLX5rOzs7Bg4cCA7\nduzIlCwAoqKiKFu2rGHugpOTEy4uLgB4enoSFxcHwNGjR2nWrJnhdSdPnqRhw4ZUr16d7777DoDb\nt2/TpEkTAgMD8fX1NbyRr1y5Ej8/P3x9ffnggw8Mx7C1teXjjz8mICCABg0acPfuXSBzD+by5cu0\nbNmSgIAAateuzZUrV7Kcn62tbbb179+/H4Bt27bRsGFDateuzWuvvUZycrLh/EJDQ6lduzZ+fn5c\nuHCBiIgI5s+fz8yZMwkKCmL//v307duXIUOGUL9+fcaOHcuJEyeoX78+/v7+dOnShfj4eODhJZVR\no0YZ6j9y5AgAS5YsYfjw4QDcuXOHV155hYCAAAICAnLck5EehpFUIVWIQy8womISVcM+c9XGPWdN\nHUqRc//+fRUQEKCqV6+u3n77bbVnzx7Dc56enio2NlYppdSRI0dUSEiIUkqpCRMmKH9/f5Wamqpi\nYmKUh4eHunXrlpo+fbr64osvlFJK6fV6lZSUpCIjI1XFihVVTEyM0mq1qnnz5uqXX35RSiml0WjU\npk2blFJKvf/+++rzzz9XSikVGhqqZsyYoZRSKjg42FA+LS1NpaSkZDkHW1tbpZTKVL9Op1NJSUkq\nOjpaNWnSxPC6KVOmqE8//dRwfnPmzFFKKfXNN9+oAQMGZKlfKaX69u2rOnbsqPR6vVJKKV9fX7V3\n716llFLjx49Xo0aNUkopFRISogYNGqSUUmrv3r2qVq1aSimlFi9erIYNG6aUUuq1115Ts2bNMsSY\nkJCQo3+nS9ejVfNB83NUtqh4nvdOSa3FmHMZO74b/yqh327n5MXbpg6nSLGxseHYsWMsWLCAcuXK\n0b17d5YuXfrU12g0Gl5++WWsrKwoU6YMzZo1IywsjODgYBYvXszEiRM5deoUtra2HDlyhGbNmlGm\nTBnMzc1544032Lt3LwAlSpSgffv2ANSuXZuIiIhM9dy/f59bt27RuXNnQ/mSJUs+Ma5H6z99+jS2\ntrYcOnSIs2fP0rBhQwIDA1m2bBnXr183vKZLly4ABAUFZapfPXb/f7du3dBoNCQkJJCQkEDjxo0B\n6NOnj+F8AHr27AlA48aNSUxMJCEhIdNxdu3axdChQwEwMzOjdOnST23rf0gPwzjSUsVcjSrOTBre\nlrcnrSf63n1Th1NgHTp0KNMbYk6YmZnRtGlTQkNDmTNnDmvXrgUeXpbV6x/ecJCa+vR5MWZmZjRu\n3Jh9+/bh5uZG3759Wb58eZbJV0opwxjIo0t4mJmZodVqjYr7cdnVD9CqVSvCw8MJDw/nr7/+YuHC\nhYbXWFk9vG3b3Nz8qfWXKpX9DOvHE8vjsttQ6t9ekx2dTmFeADanKiykpQStG1SnWys/3p70C+kZ\nsu7U4xYvXkynTp34+++/c/yaixcvcunSJcPj8PBwPD09gYfX+I8ePQpgSCLw8A1vw4YNpKWlERsb\ny+7du6lbty7Xr1+nXLlyDBgwgAEDBhAeHk5wcDB79uwhNjYWnU7HTz/9RNOmTZ8ak1IKpRS2tra4\nu7sbJsympaXx4MGDJ74uu/rr16/PgQMHDG2SnJyc6XyzY2dnR1JSUrbP2dvb4+joaBgfWb58ueFW\nUKUUq1atAmD//v04ODhkuW25RYsWzJs3D3g4gP7PxN9/o9NLD8MYciO+AGBEz0ac+TuKL5fu5qMB\nLUwdToGg1Wr5z3/+w5YtW9izZ0+m5Wv+zf379xk+fDjx8fFYWFjg5eXFggULAJgwYQL9+/endOnS\nmXbV02g0+Pn50axZM2JiYhg/fjwVKlRg2bJlfPnll1haWmJnZ8eyZcuoUKECU6ZMoVmzZiil6NCh\ng+GOoUfvttJoNJmO/8/Py5cvZ/DgwYwfPx5LS0vWrFljSGiPvhYeXu6ZPn16pvrLli3LkiVL6Nmz\np+Hury+++AIvL68sx/jnOB07dqRr1678+uuvfP3111liXbp0KUOGDCElJYWqVauyePFiQxlra2uC\ngoLQarV8//33WY49a9YsBg0axKJFizA3N+fbb7/N0VLsxWEeRm6S5c2FQXzSAzoMX8ynb7emeXA1\nU4djUgkJCfTs2ZOMjAxWr16No6OsYGoqzZo1Y8aMGQQFBeX6sY/8dYOpi3ezZnrRXFwxO4VqLSlR\ncDnYlWTmmI6MnfUbUTHZXzooLlavXk3lypXZsmWLJIsiTKeTeRjGkB6GyGL2ygP8efIaP3zRo9j+\nZ/rnb0smNRZtB05E8M3qg6yY1NPUoeQb6WGIXPX2aw9n4n6zOveWcShsHr0+LoounU5fbD8UPQtp\nKZGFubkZX43pyLJNxwk7c+PfXyCeS0hICMeOHXtqmQ0bNnDu3DnD4wkTJrBjx468Dq3I0+r0mMug\nd45JwhDZci5jx9SR7Rg9YxPJD9JNHU6emjdvHhcvXjRZ/Tnpzaxfv56zZ88aHk+cOJEWLeRutuel\n00sPwxjSUuKJmgdXo4FfRaYt3WPqUPKEUorx48fz1VdfYW2duzsRJicn0759ewICAvD19WX16tXs\n2LGDoKAg/Pz86N+/P+npWRPxP+s3AaxZs4Z+/fpx8OBBNm7cyJgxYwgKCuLKlSv07dvXMIfjScfN\nbk0nkZlOr7CQiXs5Ji0lnurjAS34/c8LRe7SlF6vZ9SoUWzatIl9+/ZRsWLFXD3+1q1bcXNz48SJ\nE5w+fZo2bdrQr18/Vq9ezalTp9BqtYaJZo96fA4FQIMGDejUqRPTp0/n+PHjVKlSxdArSU1NfeJx\nNRoN5cqV49ixYwwdOpTp06fn6jkWBTKGYRxpKfFU9nbWfPZ2az6YtYXUtAxTh5MrdDodgwcP5siR\nI+zcuZPy5cvneh1+fn5s376dDz74gP379xMREUHlypWpVu3h/JbH10rKicfvbFFKceHChace90lr\nOomHZAwm0INjAAAgAElEQVTDOJIwxL9qVb86L1SrwMwV+00dSq74448/uHz5Mtu2bcPBwSFP6vDy\n8iI8PBxfX18+/vjjLPvWP+m2xkd7GI8v15HdOMfjv3t0TSnI+ZpOxZVOFh80irSUyJHQIa1Yt+MM\nJy7cMnUoz61NmzZs374903hBbrt9+zbW1ta88cYb/Oc//+HgwYNcu3bNsPbSo2slPcrZ2Znz58+j\n1+tZv3694c3fzs4uy/pIGo0Gb29vIiIiMh3339aUEv9PKxP3jCJrSYkcKWNfivGDWjB21hZ+ndUX\nK8vC/adjYZG38Z8+fZoxY8ZgZmZGiRIlmDdvHvHx8XTr1g2tVktwcDBDhgzJ8ropU6bQoUMHypUr\nR506dQybEvXo0YOBAwcye/Zsfv75Z0N5KysrFi9enO1xn7SmlPh/erlLyigy01vkmFKKgZ+upXZN\nN4Z2y7rNphCFzQ+bj3Pu6l2+GNbW1KHkG5npLfKFRqPhk0EtWLgujNsxOVs+2tT0ej03bhStO7xE\n7pExDONISwmjVHJx5I2XApm8aJepQ/lXSimGDx/OyJEjTR2KKKC0ehnDMIa0lDDa290acPxcJIdP\nG7cDXX5SSvHBBx8QFhbGkiVLTB2OKKCkh2EcaSlhtJLWlnw4oDmh325Hq9ObOpxsTZ48mS1btrB1\n69Yc7+8sih+dXi9btBpBWko8k3aNvHGyL8WKzcdNHUoW3377LYsWLWLbtm2UKVPG1OGIAkwm7hlH\nEoZ4JhqNhgmDW/L1ygPExCebOpxMrKys2LZtGy4uLqYORRRwsoGScaSlxDOrXqkcrzSvxfRlxi1x\nkdf69etH1apVTR2GKAR0ehnDMIa0lHguw3s04o9Dl7h0PcbUoQhhNFl80DjSUuK52NtZM7hrfaYt\nKZpLoIuiTcYwjCMJQzy33h2COHflDkf+yv8JcnFxcZw6dSrf6xVFg/QwjCMtJZ6bVQkLRvdqzOTv\nd+Xrci2pqal07tyZVatW5VudomjR6pWMYRhBWkrkis4hL/AgTcu2g5fypT69Xk+/fv2oUKECn332\nWb7UKYoevfQwjCItJXKFubkZH/QLYdrS3fkymW/8+PFERESwbNkyzGTilXhGWr1etmg1grSUyDVN\ngipToYwdq7edzNN6li5dysqVK9mwYQMlS5bM07pE0abT6TEzl0HvnJKEIXKNRqNhbL8QZv14gJTU\n9Dyrx9ramk2bNuXJ1qqieNHqlPQwjGCSltq6dSs+Pj54eXkxderUbMvs3r2bwMBAatWqle3OZKJg\n8vNyoU5NN37YHJ5ndXTv3p0aNWrk2fFF8SEbKBkn31tKp9MxbNgwtm7dytmzZ1m5ciXnzp3LVCY+\nPp533nmHjRs3cubMGdasWZPfYYrnMKLniyxcd5jkB3nXyxAiN2hltVqj5HtLhYWFUa1aNTw9PbG0\ntKRHjx5s2LAhU5kff/yRV199FXd3dwDKli2b32GK5+DtWY56vhVZvqngLUwoxKN0Oj1mMnEvx/J9\nY+bIyEg8PDwMj93d3Tl8+HCmMpcuXSIjI4NmzZqRlJTEyJEj6dWrV5ZjhYaGGn4OCQmRS1cFyMjX\nX+T1cT/yZvtAbEtZPdexbt68afjwIERu0uqK/jyM3bt3s3v37lw5Vr4njJxsRJ+RkcHx48fZsWMH\nKSkpNGjQgPr16+Pl5ZWp3KMJQxQsXhXL0tDfk6Ubj/FO94bPfJwVK1YwadIkTp06hbm5eS5GKETx\nGMN4/MP0xIkTn/lY+d5Sbm5umfZYvnHjRpZPjx4eHrRu3ZqSJUtSpkwZmjRpwsmTeXurpsh9I3o2\n4vtfjpCUkvZMrz969Cjvvvsuq1atkmQh8oSMYRgn31uqTp06XLp0iYiICNLT01m1ahWdOnXKVKZz\n587s378fnU5HSkoKhw8fpmbNmvkdqnhOVT3K0KR2FZb8etTo1965c4cuXbowf/58atWqlQfRCfFw\neXMZw8i5fE8YFhYWzJkzhzZt2lCzZk3DLZLz589n/vz5APj4+NC2bVv8/PyoV68eAwcOlIRRSA3v\n2YglG46SmJya49ekp6fTrVs3+vXrxyuvvJKH0YnirjiMYeQmjcrP1eJykUajydeF7sSzGz1jE54u\nDox4/cUclQ8PD2f69OksX75clv0QearbmB8Y06cpwbU8/r1wEfE8753yv1HkuWE9GrJ04zHu53As\nIzAwkBUrVkiyEHlOdtwzjrSUyHNV3Jxo4F+JH387YepQhMhE9sMwjrSUyBfvvNaARb8cIS1da+pQ\nhDCQHfeMIwlD5IsaVZzxrVaBn7fL7nii4NDplPQwjCAtJfLN26814Ns1h8nQ6jL9/tNPP82yPIwQ\n+UHGMIwjLSXyTVANNypWsOfX3WcNv9u4cSMLFy6kfv36JoxMFFcyhmEcaSmRr97p3pBvfj6ITqfn\nypUrDBgwgNWrV+Ps7Gzq0EQxpNXLGIYxnithXL9+nSNHjnD9+vXcikcUcQ39K2FXyoqNe87w6quv\n8vHHH9OgQQNThyWKKRnDMM4zLz44f/580tLSsLW1Zd++fZiZmTFq1KjcjE0UQRqNhmHdG/Le1JXU\n8PZh2LBhpg5JFGNanezpbYxnThhVq1alZcuWhse7du3KlYBE0dc8uBoV3V0JrtskR6sXC5FXisNq\ntbnpmRNG6dKl+c9//kNKSgr29va89NJLuRmXKMLMzDR8F9qdzu8upW6tSjQOqmzqkEQxJavVGsfo\nhDFjxgy2bNlCVFQUnTp14quvvqJEiRJ5EZsowpzL2DHr/U4Mm7yB9TN7417e3tQhiWJIp1My6G0E\no1Ort7c3O3bs4MyZM7Rs2ZLPP/88L+ISxUC9WhUZ3LUe70xaLzPAhUlo5bZaoxjdUlFRUWzZsoXk\n5GRatGhB3bp18yIuUUz0f7ku7s4OhM7fbupQRDEkE/eMY/QlqRs3bhAfH8/ixYuJjY1Fq9WSkJBA\nZGQkY8eOzYsYRRGm0WiYOrIdr4xexqptJ+ne2t/UIYliRCbuGcfolurcuTMNGjTg559/ZufOnSxe\nvBilFFu2bMmL+EQxYFvKim8/6sK0xbs5fSnK1OGIYkIphU4vYxjG+NcNlC5cuICZmRleXl5PPdDt\n27dxcXHJ1eCeRjZQKnq27D/P5O938etXfXEsXdLU4YgiTqvT4/3yl/y9sXhdGcnTDZSqVq3KtWvX\nmDt3LvPmzePo0ez3Z87PZCGKppde9KFdI29GffkrOp3e1OGIIk6nl0l7xjJ6i9awsDCOHTuGXq/H\n29ubkJAQLCyeeTrHM5MeRtGk1enp9dFPBNfy4N03G5s6HFGEpaSmU+f12Zxd956pQ8lXz/PeafQ7\nfXBwMMHBwQCcP3+eRYsWkZ6ejpubG23atMHGxuaZAhECwMLcjK/HdqbzqCX4V3eheXA1U4ckiijZ\nPMl4z9U18PHxwcfHB4Bbt26xadMmunfvniuBieKrnKMNs8d2ZsgX61g97U0quzmZOiRRBMnCg8Yz\nurVmzJhBixYteOGFF/jwww/JyMgAwNXVVZKFyDW1a7rz7puNGfTZWpJS0kwdjiiCZFkQ4z3XTO8W\nLVrw2Wef5UVcQvB6u0Dq+1bk3S83yiC4yHWy8KDxZKa3KNA+GdSSpJQ0Zv6wz9ShiCJGehjGk5ne\nokArYWnONx++wsvvLsWncnk6NKlh6pBEEaHTKcxk0NsoMtNbFHhl7Esx/+MuTJi3jTOXZSa4yB1a\nWUfKaEb3MIKCgjI9rlq1apbNlITIbTWrOPPp220Y8vk61s/sQzlHuX1bPB+9Xo+5TNwzSq61lsz0\nFnmtfWMfXmlRi7cnrSc9Q2fqcEQhJ2MYxpPWEoXKu280xsGuJOO/+V1m+ovnImMYxpOEIQoVMzMN\nM//TgVOXoliw9rCpwxGFmPQwjCetJQod21JWLArtytKNx9i877ypwxGFlMzDMJ60liiUXMqWZuH4\nroz/5neOn4s0dTiiEJIehvFM0lpbt27Fx8cHLy8vpk6d+sRyR44cwcLCgnXr1uVjdKKweKGqM1++\n254hX6zj+u14U4cjChkZwzBevicMnU7HsGHD2Lp1K2fPnmXlypWcO3cu23Jjx46lbdu2Mrgpnqh5\ncDWG9WhIv9DVJCSlmjocUYhID8N4+b6RRVhYGNWqVcPT0xOAHj16sGHDBmrUyDyDd/bs2XTt2pUj\nR4488VihoaGGn0NCQggJCcmDiEVB17tDba7fjmfIF+tY+ll3SliamzokUQjoiskYxu7du9m9e3eu\nHCvfE0ZkZCQeHh6Gx+7u7hw+fDhLmQ0bNrBz506OHDmCRpN9t/HRhCGKt3FvNWPopPWMm/0b099t\n/8S/GSH+oSsmPYzHP0xPnDjxmY+V762Vk//Io0aNYsqUKYadoeSSlPg35uZmfDWmI5evxzBzxX5T\nhyMKAZ3M9DZavvcw3NzcuHHjhuHxjRs3cHd3z1Tm2LFj9OjRA4CYmBh+++03LC0t6dSpU77GKgqX\nUtYlWBTajW5jfqCcow292gf9+4tEsaXVKdlxz0j5njDq1KnDpUuXiIiIwNXVlVWrVrFy5cpMZa5c\nuWL4uV+/fnTs2FGShciRsg42LPusO6+9vwKn0iVp31hWtxXZ0+mKxxhGbsr3hGFhYcGcOXNo06YN\nOp2O/v37U6NGDebPnw/A4MGD8zskUcR4VHDg+9Cu9Pp4FQ52JWkU4GnqkEQBpJPVao2mUYV0gOCf\n8Q0hnuTwmeu8M+kXFk98DV+vCqYORxQw63acZv+Ja/z3vQ6mDiVfPc97p6RXUWTVq1WRScPb0n/i\nz1y+EWPqcEQBI2MYxpOEIYq01g2qM+6tZvT6aBURt+6ZOhxRgMgYhvGktUSR90rzWozo2Yg3P/qJ\nyLsJpg5HFBAyhmE8aS1RLPRsF8Bbnevw5kc/cTfuvqnDEQWA9DCMJ60lio23Xq7Lqy19efOjn4hN\nSDF1OMLEtHqFhYxhGEUShihWhnVvSJuG1Xlj3Epi4pNNHY4wIelhGE9aSxQ7o99sTNtG3rw+biXR\n9+TyVHGl1cnSIMaS1hLFjkajYdQbL9KhsQ89PviRO7FJpg5JmIAMehtPWksUWyNef5EuzWvR84Mf\niYqRpFHc6HRKLkkZSVpLFGvvdG9Ij7YBdB+7QnbtK2Z0Or1M3DOSJAxR7A16tR4DuwTTfewKLkRE\nmzockU+0MuhtNGktIYA32wcxrn8z3vzoJ8LPR5o6HJEPZAzDeNJaQvxPp6Y1mTbqJQZ+upZ94VdN\nHY7IY3JbrfGktYR4RLO6VZn30SuMnr6RjXvOmjockYe0ell80FiSMIR4TN0XPFj+RQ+mLt7NvJ8P\nyjL6RZT0MIwnrSVENnw8y7N2Ri827T3Hh3O2otXpTR2SyGU6vR4LmbhnFGktIZ7AuYwdq6a+QVRM\nEgMmruF+SpqpQxK5SHoYxpPWEuIpbEtZsXB8V9zKl6bbmB+4ESVzNYoKrU5hYS5jGMaQhCHEv7Aw\nN+Pzd9rQo40/Xd5bzp8nr5k6JJELdDo9ZnJJyijSWkLkgEajoU+nOsx6vxMjp/3K0o3HZDC8kNPq\nZB6GsaS1hDBCQ/9KrJ3ei5+2nmDc17+Rlq41dUjiGen1MoZhLGktIYxU0cWBNdN7kZSSRtcxy2UN\nqkJKxjCMJwlDiGdgU7IEcz54ma4tfOny3jJ+//OiqUMSRtLpZQzDWBamDkCIwuqfcQ1/b1eGT93A\n4TPX+aBfM0pYmps6NJEDMoZhPGktIZ5TgLcrG2f15XpUPN3GLOdKZJypQxI5oNfLjnvGktYSIhc4\n2JVk4Sev8mpLX7r+Zzkrfzshd1EVcA/HMOQt0BjSWkLkEo1GQ+8OtVk97Q1W/BbOwM/WEhOfbOqw\nxBM8nOktg97GkIQhRC6r5lGWdTN64+VRlvbDF/PbgQumDklkQ6uTS1LGktYSIg+UsDRnbL8Q5o57\nmS+X7uGdyeuJvie9jYJE5mEYT1pLiDxUp6Y7W2b3o5KLIy8NW8S6nWdkbKOAkLukjCetJUQes7ay\n5P2+IXwf+hoL14XRZ/xqrsqdVCan0ykZwzCSSRLG1q1b8fHxwcvLi6lTp2Z5fsWKFfj7++Pn50ej\nRo04deqUCaIUInf5elXg16/60DiwMq/+Zzkzlu/lQWqGqcMqtmQMw3j53lo6nY5hw4axdetWzp49\ny8qVKzl37lymMlWqVGHv3r2cOnWKTz75hEGDBuV3mEL8q2e5tGRpYc7ALsFsmf0WEbfu0frt7/jj\n0CW5TGUCOhnDMFq+t1ZYWBjVqlXD09MTS0tLevTowYYNGzKVadCgAfb29gDUq1ePmzdv5neYQjzV\nH3/8gYuLC++99x5//fWX0a+vUNaO2WM7M3l4O6Ys3k2vj3/i3JU7eRCpeBLZcc94+b40SGRkJB4e\nHobH7u7uHD58+InlFy1axEsvvZTtc6GhoYafQ0JCCAkJya0whXiqZs2asXbtWhYtWkSTJk2wt7dn\nwIABfPjhh0Yd58VAT36b+xY//X6S3p+splndKrzXqwnOZezyKHLxj+IyhrF79252796dK8fK94Sh\n0eT8H2jXrl18//33HDhwINvnH00YOfHFF1+wcuVKzM3NMTMzY/78+QQHBxt1DFM4duwYy5YtY9as\nWTkqP3DgQEaPHk2NGjVypf5r167x559/0rNnz1w53q1btxg5ciQ///xzrhwvN+l0OqKjo7l16xa3\nb9/O9P3Rn6Ojo3FwcMDFxYWyZcty5coVbt++/Ux1WlqY06t9EJ1DavLtz4do+84i3nwpiIFdgilt\na53LZyj+UVzGMB7/MD1x4sRnPla+Jww3Nzdu3LhheHzjxg3c3d2zlDt16hQDBw5k69atODo6Pne9\nBw8eZPPmzYSHh2NpaUlcXBxpaXm3R/M/16SNSZBPUrt2bWrXrp3j8gsXLnzuOh919epVfvzxx2wT\nhlarxcLCuD8jV1fXfE8W/ySC7N78H/1ddHQ0jo6OuLi44Orqavju5+dH27ZtDY+dnZ05deoUgwYN\nomLFimzZsoWqVas+V4ylbax5v28Ir78UyJyVB2g2cD79Xq5L3461sS1llUstIf4he3o/A5XPMjIy\nVJUqVdTVq1dVWlqa8vf3V2fPns1U5tq1a6pq1arq4MGDTzyOsaGvW7dOdezYMdvnKlWqpGJjY5VS\nSh05ckSFhIQopZSaMGGCevPNN1WDBg2Ul5eXWrhwoeE106ZNU3Xr1lV+fn5qwoQJSimlrl69qqpX\nr6569+6tXnjhBXXt2jXVp08fVatWLeXr66tmzpyplFIqPDxc1atXT/n5+alXXnlF3bt3TymlVNOm\nTdXYsWNVcHCwql69utq3b59SSqldu3apDh06KKWUSkpKUn379lW+vr7Kz89PrV27Nsv5NG3aVB07\ndkzpdLps6798+bJq27atql27tmrcuLE6f/68UkqpPn36qBEjRqiGDRuqKlWqqDVr1iillKpXr56y\nt7dXAQEBaubMmWrJkiWqY8eOqnnz5iokJETFxcWpzp07Kz8/P1W/fn116tSpp7bf1atXVa1atZRS\nSmm1WvXee++pWrVqKT8/PzV79myj/l21Wq26ffu2On78uNq0aZNasGCBmjhxoho8eLDq1KmTqlOn\njnJ1dVWWlpaqfPnyyt/fX7Vr10699dZb6uOPP1bffPONWr9+vTp8+LC6fv26SktLy1G9er1eNWnS\nRC1btkzp9XqjYs6pKzdj1agvf1V1Xp+l5q85pO6n5Cw2kTMB3WequIQUU4eR757nbT/fexgWFhbM\nmTOHNm3aoNPp6N+/PzVq1GD+/PkADB48mE8//ZR79+4xdOhQACwtLQkLC3uuelu3bs2nn36Kt7c3\nLVu2pHv37jRp0gR4ei/gzJkzHDp0iPv37xMYGEj79u05ffo0ly9fJiwsDL1eT+fOndm3bx8eHh5c\nvnyZ5cuXExwczLFjx7h16xanT58GIDExEYDevXszd+5cGjduzIQJE5g4cSIzZ85Eo9Gg0+k4fPgw\nv/32GxMnTmT79u2Z4vnss89wdHQ03GocH591855/zic8PDzb+gcNGsT8+fOpVq0ahw8f5u2332bH\njh0AREVFceDAAc6dO0enTp149dVXmTp1KtOnT2fjxo0ALFmyhPDwcE6fPo2DgwPDhw+ndu3a/PLL\nL+zatYvevXsTHh6ebft16NAhU6wLFizg+vXrnDx5EjMzM+7du5fp+Rs3bnDq1Kkc9wj+6RX4+/vT\ntm1bw++cnZ2xtLR86t+IMTQaDXv27Mm142WnspsTM//TkYvXopn14wEWrD1M7w5B9O5YGwe7knla\nd3FQXMYwcpNJ9sNo164d7dq1y/S7wYMHG37+7rvv+O6773K1ThsbG44dO8a+ffvYtWsX3bt3Z8qU\nKfTp0+eJr9FoNHTu3BkrKyusrKxo1qwZYWFh7Nu3j23bthEYGAhAcnIyly9fxsPDg0qVKhnGRapW\nrcqVK1cYMWIE7du3p3Xr1iQkJJCQkEDjxo0B6NOnD926dTPU2aVLFwCCgoKIiIjIEtOOHTtYtWqV\n4bGDg8MT48+u/vv373Pw4MFMdaanpxvO9+WXXwagRo0a3Lnz8K4d9dgtnxqNhlatWhnqPnDgAOvW\nrQMeDgbHxsaSlJSUbfsdPnwYf3//TOczdOhQw0Y2j19+nD17NmfOnMm3RFAQVa9UjrnjXubKzVjm\nrzlMs4Hz6dbKj7c616VCWRkcf1Yy09t4xWoDJTMzM5o2bUrTpk3x9fVl6dKl9OnTBwsLC/R6PQCp\nqalPPcY/n97HjRuXZX5IREQENjY2hscODg6cOnWKrVu38u2337J69WpmzpyZ6TWPvxlbWT28Vm1u\nbo5Wm/1+0Y+/5kkcHBw4efIkv//+u6H+r776CgcHB0MP4HElSpTIUT2PnqcxMWW3w9nTXjtt2rQc\nHbc4qOJehqmjXmLE3UZ8t/4Ibd9ZRJOgyvTtVIdAH9dcGS8rTnTFZNA7NxWb1rp48SKXLl0yPA4P\nD8fT0xMAT09Pjh49CsDatWsNZZRSbNiwgbS0NGJjY9m9ezfBwcG0adOG77//nuTkh4vJRUZGEh0d\nnaXO2NhYtFotXbp04bPPPiM8PJzSpUvj6OjI/v37AVi+fLlRtwO3atWKuXPnGh5nd0nqn9hjY2PR\n6XSZ6rezs6Ny5cqsWbPGUO7fZtLb2dmRlJSU6diPaty4MStWrAAe3sJXrlw57Ozssm2/unXrZjmf\n+fPno9PpALJckhJZuZW3Z8Lgluz9fgiBPq6MnrGRl99dyrqdZ0jLyP5DhshKJu4Zr9j0MO7fv8/w\n4cOJj4/HwsICLy8vFixYAMCECRPo378/pUuXJiQkxPBJTaPR4OfnR7NmzYiJiWH8+PFUqFCBChUq\ncO7cORo0aAA8fEP94Ycf0Gg0mT7lRUZG0q9fP0PvZcqUKQAsXbqUIUOGkJKSQtWqVVm8eHG2MT96\nrH9+/vjjj3nnnXfw9fXF3Nyc0NBQw2Wkx1/7pPpXrFjB0KFD+fzzz8nIyKBnz574+fk9sU5/f3/M\nzc0JCAigb9++ODo6ZioXGhrKW2+9hb+/PzY2NixduvSp7RcREWF4/YABA7h48SJ+fn5YWloyaNAg\n3n777X/75xQ8vKuqX+e69O5Qm93HrrDk16NMXrSLnu0CeK21H+7l7U0dYoGllEKnV5ibSa/MGBqV\n02sJBYxGo8nz5RQmTpyIra0t7733Xp7WU1RJ++W/yzdiWL7pOL/uOUutqhXo1tqP1vW9sLYq2uM8\nxtLq9Pi8/CWXN441dSj57nneO6U/9i/kuvDzkfbLX9U8yjJxaGsOLRvGa639WPPHaRr0mcvHc3/n\n5MXbsmbV/8j4xbORHoYQRVzk3QTW7TzDmu2nsSphQYfGPrRvXIOqHmVMHZrJJD9Ip+4bszm7rvj1\nfp/nvVMShhDFhFKK4+ci2bz/PFv2n8fRriQv/S95VHFzMnV4+SoxOZUX+87j1M/vmjqUfCeXpHJZ\nSEgIx44de2qZDRs2ZFqWfcKECYbJb0IURBqNhto13Rk/qCV/LnmHT99uTWx8Cj3GruClYd8ze+UB\n/vr7TrH4ICaT9p5NsblLyhiP3+2UnfXr19OxY0fDAn/Ps6CXEPnNzExD3Rc8qPuCB58MbMHRszfZ\nfugSw6b8QmqalmZ1q9I8uCqN/D0paV30BsyLy8KDua3YtFhycjLt27cnICAAX19fVq9ezY4dOwgK\nCsLPz4/+/fsbZjw/ytbW1vDzmjVr6NevHwcPHmTjxo2MGTOGoKAgrly5Qt++fQ1zOJ50XE9PT0JD\nQ6lduzZ+fn5cuHAhf05eiKcwNzejnm9FPh7Ygl0LB7Nick+quDvx/S9HCH5zNm+F/sziDUe4EBFd\nZHofsvDgsyk2LbZ161bc3Nw4ceIEp0+fpk2bNvTr14/Vq1dz6tQptFot8+bNy/K67OYlNGjQgE6d\nOjF9+nSOHz9OlSpVDL2S1NTUJx5Xo9FQrlw5jh07xtChQ5k+fXr+nLwQRqji5sSAV4L5cfLrHFjy\nNq+2qMXF6zEM+mwtwW/OYeS0X1m97RQ37yaYOtRnptPLsiDPoti0mJ+fH9u3b+eDDz5g//79RERE\nULlyZapVqwY8XNNp7969Rh3z8U9bSikuXLjw1OP+21pRQhQkpW2tad+4BpOHt2PPoiGs/29vGgZU\nYn/4VV55dylN+3/Le//dxE9bT3Dpegx6feHogcgYxrMpNmMYXl5ehIeHs3nzZj7++GOaN2+e6fkn\ndbUf7WE8ePDgic896XdKqUy/y8laUUIUVO7O9nRv7U/31v4opbh0PYZjZyMJ++sm834+RGJyKrVr\nuFOnpjuBNVx5oYpzgdzLQyvbsz6TYpMwbt++jaOjI2+88Qb29vbMnTuXa9eu8ffff1O1atUnrunk\n7OzM+fPnqV69OuvXrzfsNW5nZ2dYLvwfGo0Gb29vIiIiMh23adOm+XGKQuQrjUZD9UrlqF6pHD3b\nBUTysLMAACAASURBVABwN+4+R8/e5OjZm0xdvJvzV6NxKWeHn5cLtapVoFY15wKRRHQ6fbYLYYqn\nKzYJ4/Tp04wZMwYzMzNKlCjBvHnziI+Pp1u3bmi1WoKDgxkyZEiW102ZMoUOHTpQrlw56tSpY1hw\nsEePHgwcOJDZs2dn2j3OysqKxYsXZ3vcx8dDZBa0KGrKO9ny0os+vPSiD/DwbqRL12M4czmK05ei\n2LT3HBcionEtVxpfrwrUrOJM9UplqV6pLBXK2OXb/wkZw3g2MnFPCJGvMrQ6Ll+P4dTlKM5fjebi\ntWguXoshLUNL9Ypl8apUFu9K5fCqWJbqlcpR1qFUrieSv/6+w/tfbWbz7Ldy9biFwfO8dxabHoYQ\nomCwtDCnRhVnalRxzvT7uIQULl2P4cK1aC5di+G3/Re4cC0ahaKyqxOero54Gr47UtnVCXs762eK\nQTZPejaSMIQQBYKTfSnq+Vaknm9Fw++UUsQlPiAiMo6IW/eIuHWPHWGXibgVR0TkPSwtzank4kjF\nCg64lS+NW3n7TN9LWZfIti4Zw3g2kjCEEAWWRqOhjH0pytiXonZN90zPKaWIiU8h4lYcN+4kEHk3\ngTN/R/H7wYtE3k3gVnQi1iUsDMnDpWxpKpSxw7mMLfcSH0gP4xlIwhBCFEoajYZyjjaUc7Sh7gse\nWZ5XShGbkELk3UQi7yYQFZNEVGwS567eJSomiQBvFwAmTZrE4sWLqVChAs7Ozoavtm3bEhwcnN+n\nVaDJoLcQolhLTk4mMjKS27dvc+fOHe7cucPdu3dp1qxZlvlaAOPHj+fHH3+kTJkyODk54ejoiKOj\nIz169KBx48ZZyt+9e5f09HRKly6Nra2tyS+FyfLmQgiRTxITE4mKiiIuLo64uDju3bvHvXv3aNiw\nIUFBQVnKT5gwge+++46kpCSSk5OxtrbG1taW6dOn06tXryzlly9fzsGDBylZsiTW1tZYW1tjZWVF\n27ZtDVspP+rMmTPcunULCwsLLCwsMDc3x9zcnCpVqlC+fPks5SVhCCFEIaDX63nw4AFJSUnY2Nhg\nZ2eXpcy+ffs4deoUDx48IC0tzfD9lVdeoWHDhlnKz5o1i02bNqHT6cjIyECn06HT6fjggw945ZVX\nspSXhCGEECJH/q+9e4+Lqsz/AP4ZZhDlDgkoFyWCFDQGvJGtrPhTIzSJ1E0t1LxsZOWldtttdysl\nL2lpprm52sXW2NQkTVIkLwneQlTATLzghbgpioiIIAwz5/eH22yIyhwu58wZPu/Xa14xzMP4fXp0\nPjznec45vIESERG1OgYGERGZhIFBREQmYWAQEZFJGBhERGQSBgYREZmEgWGmUlNT5S6hVVly/yy5\nbwD715YxMMyUpf+lteT+WXLfAPavLZMlMFJSUtC9e3cEBARg0aJFd20zY8YMBAQEQKvVIisrS+IK\niYjoTpIHhl6vxyuvvIKUlBTk5ORg3bp1OHnyZL02ycnJOHv2LHJzc7F69WpMmzZN6jKJiOhOgsQO\nHjwoREZGGp+/++67wrvvvluvTVxcnLB+/Xrj827dugmXLl2q1wYAH3zwwQcfTXg0leT3wygqKoKP\nz/+uXe/t7Y1Dhw412qawsBAeHv+7paPA60gREUlK8kNSpt7M/c5AaOmbwBMRkTiSB4aXlxcKCgqM\nzwsKCuDt7X3fNoWFhfDy8pKsRiIiakjywOjTpw9yc3ORl5eH2tpabNiwAdHR0fXaREdHY+3atQCA\n9PR0ODs71zscRURE0pN8DUOj0WDFihWIjIyEXq/HlClTEBgYiFWrVgEA4uLiMGzYMCQnJ8Pf3x92\ndnZYs2aN1GUSEdGdmrxcLpHt27cL3bp1E/z9/YWFCxfetc306dMFf39/ITg4WMjMzJS4wuZprH97\n9uwRHB0dhZCQECEkJESYO3euDFU2zaRJkwR3d3ehZ8+e92yj5LFrrH9KHrv8/HwhIiJCCAoKEnr0\n6CEsW7bsru2UOn6m9E/J41ddXS3069dP0Gq1QmBgoPDGG2/ctZ3Y8TPrwKirqxMeeugh4cKFC0Jt\nba2g1WqFnJycem22bdsmREVFCYIgCOnp6UJYWJgcpTaJKf3bs2ePMGLECJkqbJ69e/cKmZmZ9/xA\nVfLYCULj/VPy2F28eFHIysoSBEEQbty4ITz88MMW9W/PlP4pefwEQRBu3rwpCIIg6HQ6ISwsTNi3\nb1+915syfmZ9aZCMjAz4+/vD19cX1tbWGDt2LLZs2VKvTVJSEiZOnAgACAsLQ3l5OUpKSuQoVzRT\n+gcodwtxeHg4XFxc7vm6kscOaLx/gHLHrlOnTggJCQEA2NvbIzAwEMXFxfXaKHn8TOkfoNzxAwBb\nW1sAQG1tLfR6PVxdXeu93pTxM+vAuNv5GEVFRY22KSwslKzG5jClfyqVCgcPHoRWq8WwYcOQk5Mj\ndZmtRsljZwpLGbu8vDxkZWUhLCys3vctZfzu1T+lj5/BYEBISAg8PDwwaNAgBAUF1Xu9KeMn+aK3\nGJZ+zoYpdfbq1QsFBQWwtbXF9u3bERMTgzNnzkhQnTSUOnamsISxq6ysxOjRo7Fs2TLY29s3eF3p\n43e//il9/KysrJCdnY3r168jMjISqampiIiIqNdG7PiZ9QzD0s/ZMKV/Dg4OxqllVFQUdDodysrK\nJK2ztSh57Eyh9LHT6XQYNWoUYmNjERMT0+B1pY9fY/1T+vj9ysnJCcOHD8eRI0fqfb8p42fWgWHp\n52yY0r+SkhLjbwEZGRkQBKHBsUilUvLYmULJYycIAqZMmYKgoCDMmjXrrm2UPH6m9E/J41daWory\n8nIAQHV1NXbu3InQ0NB6bZoyfmZ9SMrSz9kwpX+JiYlYuXIlNBoNbG1tsX79epmrNt24ceOQlpaG\n0tJS+Pj4ID4+HjqdDoDyxw5ovH9KHrsDBw4gISEBwcHBxg+aBQsWID8/H4Dyx8+U/il5/C5evIiJ\nEyfCYDDAYDBg/PjxGDx4cLM/O1WCkrcBEBGRZMz6kBQREZkPBgYREZmEgUFERCZhYBARkUkYGERE\nZBIGBhERmYSBQUREJmFgEBGRSRgYRERkEgYGERGZhIFBREQmYWAQEZFJGBhEd3j++efx1ltvyV1G\nA3l5ebCysoLBYJC7FGqjGBikOL6+vrCxscHVq1frfT80NBRWVlbGS1Q3lUqlUtyd44ikwMAgxVGp\nVPDz88O6deuM3zt+/Diqq6v5QW+iuro6uUsgBWJgkCLFxsYa7xYGAP/+978xYcKEevco3rZtG0JD\nQ+Hk5IQuXbogPj6+3nvs378fjz32GFxcXNClS5d671dWVoYnn3wSjo6OePTRR3H+/Pm71vHrYaK1\na9eia9eucHNzw4IFC4yv33l4KzU1FT4+Psbnvr6+WLx4MYKDg+Hg4IApU6agpKQEUVFRcHJywtCh\nQ413TvvVZ599Bi8vL3h6emLJkiXG7wuCgIULF8Lf3x8dO3bEmDFjcO3atXp1fv755+jatSuGDBli\n0v9not9iYJAiPfroo6ioqMCpU6eg1+uxYcMGxMbG1mtjb2+PhIQEXL9+Hdu2bcPKlSuxZcsWAMAv\nv/yCYcOGYebMmSgtLUV2dja0Wi2A2x+869evx5w5c3Dt2jX4+/vjH//4x33rOXDgAM6cOYPdu3fj\nnXfewenTpwE0fnhLpVJh06ZN2L17N06fPo2tW7ciKioKCxcuxOXLl2EwGLB8+fJ6P5OamoqzZ89i\nx44dWLRoEXbv3g0AWL58OZKSkrB3715cvHgRLi4uePnll+v97N69e3Hq1Cl8//33JvxfJqqPgUGK\nNX78eKxduxY7d+5EUFBQgxvYDxw4ED169AAAPPLIIxg7dizS0tIAAF999RWGDh2KMWPGQK1Ww9XV\n1RgYKpUKI0eORJ8+faBWq/Hcc88hOzv7vrXMnj0bNjY2CA4OhlarxbFjx4yvNXZTy+nTp8PNzQ2e\nnp4IDw9H//79odVqYWNjg6effhpZWVkN/qwOHTqgZ8+emDRpkvHQ3L/+9S/MmzcPnp6esLa2xuzZ\ns5GYmFhvkXzOnDno0KEDbGxs7lsT0d2Y9T29ie5FpVJh/PjxCA8Px4ULFxocjgKAQ4cO4Y033sCJ\nEydQW1uLmpoaPPPMMwCAgoIC+Pn53fP9PTw8jF936NABlZWV962nU6dOxq9tbW0bbX+/P+u3z9u3\nb9/gvX57SKtLly44fvw4gNuzpqeffhpWVv/7PVCj0aCkpOSuP0skFmcYpFhdunSBn58ftm/fjpEj\nRzZ4/dlnn0VMTAwKCwtRXl6OF1980RgqXbp0wblz51q9Rjs7O1RVVRmfX7p0qdGfaWxG8ttdYPn5\n+caZVZcuXZCSkoJr164ZH1VVVejcubOxPTcFUHMwMEjRPvvsM/zwww/o0KFDg9cqKyvh4uKCdu3a\nISMjA1999ZXxtWeffRa7du3Cxo0bUVdXh6tXrxoPIzX2gS1GSEgIkpOTce3aNVy6dAkffvhhs99z\n3rx5qK6uxokTJ/DFF19gzJgxAIAXX3wRf//7342BcuXKFSQlJTX7zyP6FQODFM3Pzw+9evUyPv/t\nb9Aff/wx3n77bTg6OmLu3LnGD1bg9m/jycnJWLJkCR544AGEhobip59+Mr7Hnb+JN7ZwfS/jx4+H\nVquFr68vnnjiCYwdO7bR3/J/+/qdtahUKgwcOBD+/v4YMmQIXn/9deOOp5kzZyI6OhqPP/44HB0d\n0b9/f2RkZJhUJ5EpVEJL/jpFREQWS/IZRkFBAQYNGoQePXqgZ8+eDbYMAre3DTo5OSE0NBShoaGY\nN2+e1GUSEdEdJN8lZW1tjaVLlyIkJASVlZXo3bs3hg4disDAwHrtBg4cyOOvRERmRPIZRqdOnRAS\nEgLg9olVgYGBKC4ubtCOR8qIiMyLrOdh5OXlISsrC2FhYfW+r1KpcPDgQWi1Wnh5eWHx4sUICgpq\n0IaIiMRr8i/kgkxu3Lgh9O7dW9i8eXOD1yoqKoSbN28KgiAIycnJQkBAQIM2MpYuidmzZ8tdQquy\n5P5Zct8Egf1TuuZ8dsqyrVan02HUqFGIjY1FTExMg9cdHBxga2sLAIiKioJOp0NZWZnUZRIR0W9I\nHhiCIGDKlCkICgrCrFmz7tqmpKTEOGXKyMiAIAhwdXWVskwiIrqD5GsYBw4cQEJCAoKDgxEaGgoA\nWLBggfHs1Li4OCQmJmLlypXQaDSwtbXF+vXrpS5TdhEREXKX0KosuX+W3DeA/WvLFHvinkql4k4q\nIiKRmvPZyUuDEBGRSRgYRERkEgYGERGZhIFBREQmYWAQEZFJGBhERGQSBgYREZmEgUFERCZhYBAR\nkUkYGETU5tTo6nCtolruMhRH1vthEBFJRVenx/6sPGzddxI703MRHNAJCfPHyV2WojAwiMhiCYKA\noyeLkJSag+T9p9DV0wXRvw/E//X1x5oth+UuT3EYGERkcXLzS7FlzwlsSctB+3YaPDWoBzYtmYAu\nnZ0BANmni1GnN8hcpfIwMIjIIlTcvIXv0k5i486fUHK1EiMGBmLVmyMR+KB7g1s6q62soDcwMMRi\nYBCRYhkMAjJ+LsDXO45hV8ZZhIf64tXYcAwI8YVafe89PRq1CnV63h5BLAYGESnOlWs38fWOY/h6\nx0/oYGONZx4Pxpt/HAxXJ1uTfl6ttoKeh6REY2AQkSIIgoDDJwqRsC0TaZnnEfW77lj+16cQHNCp\nwSGnxmjUVlzDaAIGBhGZtZvVtfh2zwkkbMtEjU6P8cNDMe/lSDjat2/ye3KG0TQMDCIyS/kXy7Em\n6Qg2//AzHn2kC/4xdTB+F9JV9GzibhgYTcPAICKz8et5E59tzsChnwsw5nEtkldMhqebY4v+ORor\nK9Rxl5RoDAwikl2d3oDvD57Gp5sPo+x6FSbH9MXi156EXYd2rfLncYbRNAwMIpJN9S0dNuw4hk83\nH4anmwNeHP0ohoT533dLbEvgonfTMDCISHIVlbfw5bZMfJF0BL0CvbDijacQ0s1Tsj9frVZBz/Mw\nRGNgEJFkSstvYs2WI/hqexYG9XkI/1kwDg93dZO8Ds4wmoaBQUStruTqDfwr8RA2//AzRvw+EEkf\nPg+fTs6y1cM1jKZhYBBRq7lyrRIrN6Zj0+6f8Yehj2DHyqlwd7WXuyzukmoiBgYRtbjS8ptY/c0h\nfL3jJ4wc3NNsguJXVlYqGAwCDAYBVlbNP6+jrWBgEFGLKbtehU82ZWBdSjaeighCyj+noFNHB7nL\nakClUkGjvn3FWisrtdzlKAYDg4iarepWLT779jDWbDmCJx7rhm0fTYKXu5PcZd3Xr+sY1hoGhqkY\nGETUZLo6PTZ8fwwfrT+IsJ4+2LRkAnw9XeQuyyS3ZxjcWisGA4OIRBMEAcn7T2Hx2r3w9nDCp2+P\nxiMBneQuS5Tb52Jw4VsMBgYRiZL+Uz7e/fwHCAIw96VIDAj1lbukJlFb8VwMsVr3/Pu7KCgowKBB\ng9CjRw/07NkTy5cvv2u7GTNmICAgAFqtFllZWRJXSUR3yiu+hrh53+D1D7dh6tP98O3SiYoNC4C3\naW0KyWcY1tbWWLp0KUJCQlBZWYnevXtj6NChCAwMNLZJTk7G2bNnkZubi0OHDmHatGlIT0+XulQi\nwu3LeKzYcBCJu47jjyP7YflfnoJNO+UfnODZ3uJJPuqdOnVCp063j3Xa29sjMDAQxcXF9QIjKSkJ\nEydOBACEhYWhvLwcJSUl8PDwqPdec+bMMX4dERGBiIiIVq+fqK2o0xuwLiUby7/ajyGPBuD7j6fC\nzcVO7rJazO1dUpa/6J2amorU1NQWeS9Zf03Iy8tDVlYWwsLC6n2/qKgIPj4+xufe3t4oLCy8b2AQ\nUcs5kJ2H+FW74OZih7VzxyDQz6PxH1IYjVrVJmYYd/4yHR8f3+T3ki0wKisrMXr0aCxbtgz29g3P\nABWE+snfEnfZIqL7K75Sgfmf/oCfci/irT8OxtBHAyz23x6vJyWeLIGh0+kwatQoxMbGIiYmpsHr\nXl5eKCgoMD4vLCyEl5eXlCUStSk1ujp8/u1hfLIpA+OH98KS14ajvY213GW1Kq5hiCf5LilBEDBl\nyhQEBQVh1qxZd20THR2NtWvXAgDS09Ph7Ozc4HAUEbWMvUfPI+rlz3E0pwibP5iAV2PDLT4sAO6S\nagrJZxgHDhxAQkICgoODERoaCgBYsGAB8vPzAQBxcXEYNmwYkpOT4e/vDzs7O6xZs0bqMoks3qXS\nG4hftRM55y/j7bghGNzPX+6SJMUZhniSB8aAAQNgMCHVV6xYIUE1RG2PXm9AwrZMLFt3ALHDe2Hp\nn0e0iRnFnbiGIZ7yN1MTkclOni/B3z5KgY21Bl+/9xz8fTrKXZJsOMMQj4FB1AZU3arFsq8O4Jtd\nx/H68wPxhyHBbf4+EG3lPIyWxMAgsnB7j57HP/75PfoEeSPl4yno6Gw5J981R1s5D6MlMTCILFRF\n5S3M+3Q3Dh7Lx7vTn0B4rwflLsmscJeUeJJvqyWi1vdDxlk88fJnsLHWIOWfkxkWd8E1DPE4wyCy\nIOU3qvHO6t04mlOIJX96Ev2Du8pdktniLinxOMMgshA708/giZc/g5N9e2z/52SGRSM4wxCPMwwi\nhSu7XoX4VbvwU+5FLP/LU+jX06fxHyLOMJqAMwwiBdvx4+1ZhburHZI/msywEIEzDPE4wyBSoMqq\nGsxdvRvpx/Px8d+fRp8gb7lLUpzbu6R4HoYYnGEQKczRnEIMn377+mrbPprEsGginochHmcYRAqh\nq9Nj+VcHsP77Y5j/SiQe7/+w3CUpmpUV1zDEYmAQKcC5gqt4dfF3eMDZFskrJsHNpeFNx0gcrmGI\nx8AgMmOCIODLbZn48D/78VpsOJ4bFmqxd8CTmlptZdKVs+l/GBhEZurKtUr85cNklF2vQuL7sfDz\nfkDukizK7TUMLnqLwUVvIjO0L/MCnpzxBXo+1AmJi8czLFoBz8MQjzMMIjOiq9NjacI+bP7hBD58\nfQTP1m5FGisr1PGQlCgMDCIzUXj5Oma+lwRHOxt8t/x5Xoa8lXGGIR4Dg8gMpBw8jTdXfI8XRoVh\n6tP92vzNjaSgUVuhtk4vdxmKwsAgklFNbR3mf/oDUo+ew6ezRyOkm6fcJbUZarUV9DU6uctQFAYG\nkUzOFVzFK4u2wM/LFVuXT4KjXXu5S2pTeB6GeAwMIhl8u+cE5q7ehT9PGIixT2h5boUMuIYhHgOD\nSEI1ujrMW70bB7LzkLBgHAIfdJe7pDZLY6VCHS8+KAoDg0gihZev45V3v0Xnjg749sOJPAQlM84w\nxGNgEEkg7eh5/PmDrXhh1KOY+nRfHoIyA1zDEI+BQdSK9HoDlq87gK93/IR//u1p3uDIjHCGIR4D\ng6iVlF2vwqzF36FWp0fSsom8wqyZYWCIx2tJEbWC7NPFiJ75BXr4eSBh/liGhRniISnxOMMgamFf\nbc/CB1/uw7sznsDQR3mTI3N1+xatDAwxGBhELaRGV4c5/9qJozlF2Ph+LB70cpW7JLoPzjDEY2AQ\ntYCSqzfw0rvfws3ZDpuWjIe9rY3cJVEj1GoV9LwfhiiSr2FMnjwZHh4eeOSRR+76empqKpycnBAa\nGorQ0FDMmzdP4gqJxMk8WYSnXv03Ivr44eO/P82wUAjOMMSTfIYxadIkTJ8+HRMmTLhnm4EDByIp\nKUnCqoiaZn1KNhav3YtFs4ZhcD9/ucshEbhLSjzJAyM8PBx5eXn3bSMInCaSeavV6RG/aicOHS/A\n1+89xzviKZDaijMMscxuDUOlUuHgwYPQarXw8vLC4sWLERQUdNe2c+bMMX4dERGBiIgIaYqkNu1y\nWSVeWrAZLo622Lx0Ahx4CEqR1GpVm9gllZqaitTU1BZ5L5Ugw6/zeXl5GDFiBI4fP97gtRs3bkCt\nVsPW1hbbt2/HzJkzcebMmQbtVCoVZyIkuWNnLmLa/E0YE6nF9LG/442OFOxoTiHmf/YDNi259+Fx\nS9Scz06zO3HPwcEBtra2AICoqCjodDqUlZXJXBUR8F1aDibP/hqz44Zi5rMDGBYKxzUM8czukFRJ\nSQnc3d2hUqmQkZEBQRDg6sr97CQfg0HA0v/sw+bdPyNh/lgE+nnIXRK1AI3aCnpe3lwUyQNj3Lhx\nSEtLQ2lpKXx8fBAfHw+d7vZtEuPi4pCYmIiVK1dCo9HA1tYW69evl7pEIqOqW7X48wfbUFJWic1L\nJ8LNxU7ukqiFcIYhnixrGC2BaxjU2oqvVOCFud+gm68bFkx/AjbWZjchp2bIzS/FSws2Y+e//ih3\nKZJqzmcn/wUQ3UX26WK8OH8TJkX3wQujwnj/CgvEGYZ4DAyiO/x6v+1Fs4ZhSFiA3OVQK9FYWaGu\nDWyrbUkMDKL/MhgELPlyL5LScvCfd8ehuy/vt23JeC0p8RgYRABuVtfi1cXf4fqNW/h26UQ84GQr\nd0nUyngtKfHM7jwMIqkVXb6O0a8nwMWhA76cP5Zh0UZwDUM8zjCoTfsp9yJemPsNJj/VF38c2Y+L\n220IZxjiMTCozdrx4xn8bfl2LJgehcjHeGe8toYzDPEYGNTmCIKAz7ccwSebDuHz+Gegfbiz3CWR\nDLhLSjwGBrUpdXoD3lm1C4d+zkfi4vHwdneSuySSCWcY4jEwqM2orKrB9EVbUKc3YOP7sXC0ay93\nSSQjrmGIx11S1CZcLK3AH/7yH3Tq6IDP5/yBYUHGqw0beAFCk4kOjJMnT7ZGHUSt5uezlzDytS/x\n9KAeWPDKE7DWqOUuicwEZxniiA6MRYsW4dy5c61RC1GL23UoFxPf2oC344bwmlDUANcxxBG9hlFV\nVYXp06fj8uXLcHd3R79+/RAWFoa+ffti//79iImJaY06iUT7IukIVm5Mx2dz/oCQbp5yl0NmiDul\nxGnS5c0zMzNRWVkJPz8/ZGdnIyMjAxkZGcjOzsalS5dao84GeHlzuhe93oC5n+zGwWO/4PM5f4C3\nB3dC0d1pxyxF2qcvwtmhg9ylSEbyy5v36tULAPDjjz+iY8eOeOeddwAAH3/8cZOKIGopN6trMfO9\nJNyq1SHx/Vg42nNxm+5NY8U1DDGata22f//+0Ol02LdvHxwcHPDSSy+1VF1Eol0qvYGp7ySi50Me\nmPtyJBe3qVFcwxBHdGBcvXoVeXl5yM/PR35+PgoKCpCfn4/z589jwIAB+PDDD1ujTqL7Onm+BFPf\n+Qaxw0Px4uhHubhNJuEuKXFEB8aDDz6IqKgohIeHw9fXF+Hh4fDx8YGbm1tr1EfUqD2Hz+HPH2zF\nOy89juHhgXKXQwpye4bBtVBTiQ6MhQsXol+/fvjll1+g0+lw/vx5XLt2Db1798bWrVsRGxvbGnUS\n3dWX2zLx0boD+OTt0egV6CV3OaQwGjV3SYkhOjB+Xafo06eP8XsVFRU4fPgwli9fzsAgSej1Brz7\n+R6kHjmPxPfHo0tnZ7lLIgVSW3ENQ4wWuZaUo6MjBg8ejA8++KAl3o7ovqpu1WLW+9+hsqoG3ywe\nDycH7oSiptGoraDnDMNkLXotqQEDBrTk2xE1cLmsEmP/+hWc7Nvji3fGMCyoWXhfb3F4tVpSjFN5\nlzF1TiLGRYXgpWf6cycUNRt3SYnDwCBF2Hv0PF5bshVvxw1B9MAgucshC8HzMMRhYJDZ+09yFpZ9\ntR//enMk+gR5y10OWRDOMMRhYJDZ0usNWPRFKnZnnMXG92PRtbOL3CWRhVFbcdFbDAYGmaWqW7V4\ndfF3qKiswabFE7i4Ta2CMwxxeMc9Mju/7oRysG2Pf8/lTihqPVzDEIeBQWblVN5ljHxtLR5/7GG8\n/+owtLPmBQSp9XCGIQ4PSZHZSDt6Hn9ashWz44ZgBHdCkQR4HoY4DAwyCwnbMrF83QGsenMkenMn\nFEmEMwxxJD8kNXnyZHh4eOCRRx65Z5sZM2YgICAAWq0WWVlZElZHUtPrDZj/6W6sSTqCje/HUX6L\n2gAAERJJREFUMixIUtwlJY7kgTFp0iSkpKTc8/Xk5GScPXsWubm5WL16NaZNmyZhdSSlqlu1mLZg\nM34+W4JNiydw2yxJjjMMcSQPjPDwcLi43PuDISkpCRMnTgQAhIWFoby8HCUlJVKVRxL57TWhuBOK\n5MJdUuKY3RpGUVERfHx8jM+9vb1RWFgIDw+PBm3nzJlj/DoiIgIRERESVEjNVVp+E6P+/CXGRmp5\nTSiSldpKZfGBkZqaitTU1BZ5L7MLDAAQhPq7Fu71gfLbwCBlMBgEvLZkK0b8PhAvj3lM7nKojVO3\ngUNSd/4yHR8f3+T3MrvzMLy8vFBQUGB8XlhYCC8v3knNUqzc+CNu1ejw2vjfy10K0e37YVh4YLQk\nswuM6OhorF27FgCQnp4OZ2fnux6OIuU5dDwfXyQdxbK/REOjNru/etQGqdVWqDPwPAxTSX5Iaty4\ncUhLS0NpaSl8fHwQHx8PnU4HAIiLi8OwYcOQnJwMf39/2NnZYc2aNVKXSK2gtPwmZr3/HRa/Nhyd\nOzrKXQ4RAM4wxJI8MNatW9domxUrVkhQCUnFYBDw6uLvMHJwTwzs7Sd3OURGaivLX8NoSTwuQK3u\no/UHUFurx6ux4XKXQlQPZxjimOUuKbIcO348gw3fH8O3Sydy3YLMTlvYJdWS+C+YWs2ZX67gb8u3\nY+U/RsLd1V7ucogauH3xQQaGqRgY1CrKb1Tjhbnf4B9T/w/ahzvLXQ7RXWnUVqjjtaRMxsCgFlen\nN2D6oi0YEhaAkYPvfZFJIrmprbiGIQYDg1rce1+kAgLwxuRBcpdCdF8atRX0PA/DZFz0phaVsC0T\nu9Jz8c2SCVzkJrPHiw+Kw8CgFrMrPRcfrT+Ir997Di6OHeQuh6hRvLy5OPwVkFpE9uli/HVZMla9\nOZL3tSDF4AxDHAYGNdsvF68hbt4mLJo1DCHdPOUuh8hkt2cYXMMwFQODmqXsehUmvf01Zoz7HYaE\nBchdDpEoaisVb9EqAgODmqyi8hYmvLUBw8MD8dywULnLIRKNaxjiMDCoSSqravD87K/Rr6cPXhvP\na0SRMnENQxwGBolWfUuHqe98g+6+7njrj4N5i1VSLM4wxGFgkCg1ujrEzd8ETzdHzHs5kmFBisYZ\nhjg8D4NMVlNbh5cXfgsHWxu8N2sYrKwYFqRsaisVZxgicIZBJqm6VYup8YmwsdZg6Z9H8Cxusghq\ntRV3SYnAf/XUqIrKW5jw5gZ4ujli+V+i0c5aLXdJRC2C52GIw8Cg+7p6vQrj/vYVHgnojHdnREHN\nmQVZEF6tVhyuYdA9FV+pwIQ3NyBqQDe8FhvOBW6yONwlJQ5/XaS7Op57CaP+9CXGRGrxp/G/Z1iQ\nReIuKXE4w6AGdqXn4q/LkjF/+hN44rFucpdD1Go4wxCHgUFGgiBgzZYjWL3pED6Pf4a3ViWLp1bz\nWlJiMDAIAFCr02Pu6l3IOFGAxMXj4e3uJHdJRK2OMwxxGBiE4isVeGXht3jAyQ4b34+Fo117uUsi\nkgR3SYnDwGjj9mfl4bUl32HyU33xwqgwnr1NbQrPwxCHgdFGGQwCPt74I77cmollf4lG/+CucpdE\nJDnukhKHgdEGFV+pwOtLt6FWp8eWpRPRqaOD3CURyYJrGOLwPIw2JiktB9Ezv0B/bVesW/gsw4La\nNF5LShzOMNqI6zdu4a2Pv0fO+ctYE/8MHgnoJHdJRLLTWDEwxOAMw8IJgoDk/acQ+fKncHGyxXfL\nnmdYEP2XWq2CnoveJuMMw4IVXr6O2R/vQEHJdXz016fQt4eP3CURmRWuYYgjywwjJSUF3bt3R0BA\nABYtWtTg9dTUVDg5OSE0NBShoaGYN2+eDFUql65Oj082ZSB65hfoFeiFrcsnMSyI7kKjUaOuTo83\nlm/Huu3ZyDlfAl2dXu6yzJZKEARJ52N6vR7dunXDrl274OXlhb59+2LdunUIDAw0tklNTcUHH3yA\npKSke76PSqWCxKWbPUEQsPvQWSz4fA+8PZzwzrTH4evpIndZRGbt5PkSHM4pxLEzF/HTmYsoulyB\nQD93aB/ujJ7+nRDo6wY/nwdgY20ZB2Sa89kp+f+BjIwM+Pv7w9fXFwAwduxYbNmypV5gADCpQ3Pm\nzDF+HRERgYiIiBasVFlOnCvB/E93o7S8Cm+/MBgDe/vxCrNEJgj080Cgn4fx+Y2qGvx89hKOnb6I\nHzLOYuXXP6Kg5Dq6dHJGd183dPvvI6BLR3i7O5n9PWJSU1ORmpraIu8l+QwjMTER33//PT755BMA\nQEJCAg4dOoSPPvrI2CYtLQ0jR46Et7c3vLy8sHjxYgQFBdUvnDMMAMCFojJ8tP4g9mddwMxnB2BM\npJa3TyVqYTW1dThbcBWn8y7jVN4VnM67gnOFZSgtvwkfDyc86OWKB71c4eftCl9PF/h4OMPd1d4s\n/y0qaoZhym+9vXr1QkFBAWxtbbF9+3bExMTgzJkzElSnHBeKyrBiw0GkHj6H56P7IH7aUDjY2shd\nFpFFsmmnQY+HPNDjIY96379Vo0PexWu4UHQNF4rKcOREITbu+AlFlytQdr0Kbq728HJzhJe7I7zc\nneD536893R3h5mwPR3sbRR0JkDwwvLy8UFBQYHxeUFAAb2/vem0cHP53MllUVBReeukllJWVwdXV\nVbI6zdWZX65g1TeHsOfwOTw/ojf2fBIHR3teLJBIDu1trNHd1x3dfd0bvFar0+NS6Q0UXbmO4ssV\nKLp8Hdmni7Ft/0kUX65AaXkVbtXq8ICTHTq62KKjsx06OtvhAefbXz/gZAtH+/ZwsLWBg93th6Ot\nDextbWS75pvkgdGnTx/k5uYiLy8Pnp6e2LBhA9atW1evTUlJCdzd3aFSqZCRkQFBENp0WAiCgP3Z\nefhs82GcOFeCCU/2QiqDgsistbNWo0tnZ3Tp7HzPNjW1dSgtv4nS8iqUlt/E1f/+t/hKBX4+ewkV\nN2tw49dHVQ0qKm+hqkYHu/bt4GBnA3vbdrCx1sCmnQY27dRoZ62BjbUaNu00aGetho21BtbWaqit\nVFCpVM0OGskDQ6PRYMWKFYiMjIRer8eUKVMQGBiIVatWAQDi4uKQmJiIlStXQqPRwNbWFuvXr5e6\nTLNQdasW3+09iS+SjsJgMGBKTD+senMkbNpZxm4NorbOpp0GXu5O8BJx/xm93oCb1bWoqKpBZVUN\namr1qNXVoaZWjxpdHWp1etTU1qFGV4ea2jro6gwwGAQIggCDoXnrvpIvercUS170PnGuBOtSsrF1\n70n07eGD8cNDEd7rQUUd6yQi89Scz04Ghpm4cq0S2/adwuYffsbV8iqMidTiD0ODeXFAImpRDAyF\nqrh5Czt/zMWWtBwcO1OMIWEBeGpgEH4X4mv2e7uJSJkYGApScvUGdh06ix0/nkHmySI8GtwVMYOC\n8H99/dGhvbXc5RGRhWNgmLE6vQHHThdjb+YFpB09j7ziaxjU9yE83v9h/L7Xg7Dr0E7uEomoDWFg\nmBG93oBTeVdw5EQBfvwpHz8e/wXe7k4I7/Ugft/LD317eMNao5a7TCJqoxgYMrpZXYvs08U4mlOI\nwzmFyD5VjE4dHdCnhzf69fDBgFBfuLnYy10mEREABoZkrly7iVMXLiPnfAlyzl/GyQuXUVBSjh5+\nHujTwxt9grzRO9AbLo4dJK2LiMhUDIwWpKvTo6DkOvKKypBXfA2/XLyG80VlOJ13BbU6PQIfdEeQ\nnzuC/DwQ6OeOhyzossdEZPkYGCYQBAFVt3S4dqMaV8oqcelqJS5fvYGSX78uu4HCkuu4eOUGPDra\nw9fTFQ96usDX0wVdO7ugm68bPN0cefIcESlamw2M+Z/uRp3eAJ1OD53eAF2dHnV1BtTq6lBZVYuK\nm7dQcfP29VcqbtbAWqOGs0N7uLvaw+MBe3i4Otz+7wO3/+vl5ghvD2e0s+aiNBFZJkVd3rwlPeBk\nC2uNGtbWalirrYxfa9RWt6/saNcejnY2xis+MgiIiJpO0TMMhZZORCSb5nx28voTRERkEgYGERGZ\nhIFBREQmYWAQEZFJGBhERGQSBgYREZmEgUFERCZhYBARkUkYGEREZBIGBhERmYSBQUREJmFgEBGR\nSRgYRERkEgYGERGZhIFBREQmYWAQEZFJGBhERGQSBgYREZmEgWGmUlNT5S6hVVly/yy5bwD715bJ\nEhgpKSno3r07AgICsGjRoru2mTFjBgICAqDVapGVlSVxhfKz9L+0ltw/S+4bwP61ZZIHhl6vxyuv\nvIKUlBTk5ORg3bp1OHnyZL02ycnJOHv2LHJzc7F69WpMmzZN6jKJiOgOkgdGRkYG/P394evrC2tr\na4wdOxZbtmyp1yYpKQkTJ04EAISFhaG8vBwlJSVSl0pERL8lSGzjxo3C1KlTjc+//PJL4ZVXXqnX\n5sknnxQOHDhgfD548GDhyJEj9doA4IMPPvjgowmPptJAYiqVyqR2tzPh3j935+tERNS6JD8k5eXl\nhYKCAuPzgoICeHt737dNYWEhvLy8JKuRiIgakjww+vTpg9zcXOTl5aG2thYbNmxAdHR0vTbR0dFY\nu3YtACA9PR3Ozs7w8PCQulQiIvoNyQ9JaTQarFixApGRkdDr9ZgyZQoCAwOxatUqAEBcXByGDRuG\n5ORk+Pv7w87ODmvWrJG6TCIiulOTVz8ksn37dqFbt26Cv7+/sHDhwru2mT59uuDv7y8EBwcLmZmZ\nElfYPI31b8+ePYKjo6MQEhIihISECHPnzpWhyqaZNGmS4O7uLvTs2fOebZQ8do31T8ljl5+fL0RE\nRAhBQUFCjx49hGXLlt21nVLHz5T+KXn8qqurhX79+glarVYIDAwU3njjjbu2Ezt+Zh0YdXV1wkMP\nPSRcuHBBqK2tFbRarZCTk1OvzbZt24SoqChBEAQhPT1dCAsLk6PUJjGlf3v27BFGjBghU4XNs3fv\nXiEzM/OeH6hKHjtBaLx/Sh67ixcvCllZWYIgCMKNGzeEhx9+2KL+7ZnSPyWPnyAIws2bNwVBEASd\nTieEhYUJ+/btq/d6U8bPrC8NYunnbJjSP0C5O8LCw8Ph4uJyz9eVPHZA4/0DlDt2nTp1QkhICADA\n3t4egYGBKC4urtdGyeNnSv8A5Y4fANja2gIAamtrodfr4erqWu/1poyfWQdGUVERfHx8jM+9vb1R\nVFTUaJvCwkLJamwOU/qnUqlw8OBBaLVaDBs2DDk5OVKX2WqUPHamsJSxy8vLQ1ZWFsLCwup931LG\n7179U/r4GQwGhISEwMPDA4MGDUJQUFC915syfpIveovRUudsmCtT6uzVqxcKCgpga2uL7du3IyYm\nBmfOnJGgOmkodexMYQljV1lZidGjR2PZsmWwt7dv8LrSx+9+/VP6+FlZWSE7OxvXr19HZGQkUlNT\nERERUa+N2PEz6xmGpZ+zYUr/HBwcjFPLqKgo6HQ6lJWVSVpna1Hy2JlC6WOn0+kwatQoxMbGIiYm\npsHrSh+/xvqn9PH7lZOTE4YPH44jR47U+35Txs+sA8PSz9kwpX8lJSXG3wIyMjIgCEKDY5FKpeSx\nM4WSx04QBEyZMgVBQUGYNWvWXdsoefxM6Z+Sx6+0tBTl5eUAgOrqauzcuROhoaH12jRl/Mz6kJSl\nn7NhSv8SExOxcuVKaDQa2NraYv369TJXbbpx48YhLS0NpaWl8PHxQXx8PHQ6HQDljx3QeP+UPHYH\nDhxAQkICgoODjR80CxYsQH5+PgDlj58p/VPy+F28eBETJ06EwWCAwWDA+PHjMXjw4GZ/dqoEJW8D\nICIiyZj1ISkiIjIfDAwiIjIJA4OIiEzCwCAiIpMwMIiIyCQMDCIiMsn/A5Epk+BeFb9QAAAAAElF\nTkSuQmCC\n"
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.cm as cm\n",
      "\n",
      "xx, MM = np.meshgrid(x, M)\n",
      "fig = plt.figure()\n",
      "ax = fig.add_subplot(111)\n",
      "im = ax.imshow(MM.transpose(), extent=(0, 3, -1, 1), cmap=cm.RdYlBu_r)\n",
      "cb = fig.colorbar(im)\n",
      "cb.set_label(\"Mach number\")\n",
      "\n",
      "import matplotlib.patches as patches\n",
      "\n",
      "patch = patches.PathPatch(nozzle.get_paths()[0], fc='none', lw=2)\n",
      "ax.add_patch(patch)\n",
      "im.set_clip_path(patch)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAADtCAYAAACriF5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0FGW++P93dXe2zp6QhKQ7ECFAFiAJAgEhAiMBA0ME\nRE3UkYuIXK/L5Q/Pzzn3zIwX71yP3uP9jjrMePTMjI6jInOdcWRkGWWUAAkhLGENJGHJvm9k37rr\n90ekNWRr0tk6+bzOqQNV9VTVpwr48OSpp55HUVVVRQghxLihGe0AhBBCDC1J7EIIMc5IYhdCiHFG\nErsQQowzktiFEGKckcQuhJhQfHx8UBTFqsXHx2e0wx0URbo7CiEmEkVR+JiZVpV9jBzsMUXqRjsA\nIYQYaRpr2yrMwxrGsJGmGCHEhKPRWLf0prCwkBUrVhAZGcns2bN5++23e5T5+OOPiYqKYu7cuSxZ\nsoTz589b9oWEhDB37lxiYmJYuHDhsNyf1NiFEBOO1TX2Xjg4OPCrX/2K6OhoGhsbufvuu4mPjyc8\nPNxSZtq0aRw5cgRPT08OHjzI008/TXp6OtDVFHT48OFhbb+XxC6EmHB02sEfO3nyZCZPngyAm5sb\n4eHhlJSUdEvsixcvtvw+NjaWoqKibucY7nZ7SexCiAmnrxr7BVMzF00tVp8nLy+PzMxMYmNj+yzz\n+9//njVr1ljWFUVh5cqVaLVatm/fzrZt26y+nrUksQshJpy+EnuURk+Ug96y/mlHdZ/naGxsZNOm\nTbz11lu4ubn1Wubbb7/lD3/4A6mpqZZtqampBAYGUllZSXx8PGFhYcTFxQ3uRvogL0+FEBOOLS9P\nATo6OnjwwQd5/PHHWb9+fa9lzp8/z7Zt29i7dy/e3t6W7YGBgQD4+fmxYcMGMjIyhvTeQBK7EGIC\nsiWxq6rK1q1biYiIYMeOHb2WKSgoYOPGjXz00UeEhoZatjc3N9PQ0ABAU1MTX331FXPmzBny+5Om\nGCHEhGNLr5jU1FQ++ugjS5dFgFdffZWCggIAtm/fziuvvEJtbS3PPPMM0NWTJiMjg7KyMjZu3AhA\nZ2cnjz32GKtWrbLtZnohX54KISYURVH42s+6L0/jK+XLUyGEsAtarTLaIQwrSexCiAnHlqYYeyCJ\nXQgx4UhiF0KIcUYSuxBCjDOS2IUQYpyRxC6EEOPMeE/sNt3ek08+SUBAQL9fTr3wwgvMmDGDqKgo\nMjMzbbmcEEIMCZ3OusVe2ZTYt2zZwsGDB/vcv3//fq5evUpubi7vvfee5SssIYQYTRrFusVe2ZTY\n4+Liug1uc7u9e/eyefNmoGtM4rq6OsrLy225pBBC2MzWQcDGumH9YaO4uJjg4GDLutFopKioiICA\ngG7lFMWO/2sUQow4Wz/zt+ekbY1hv73b/wD6SuKqqtrN8vLLL496DBLz2FvsLV57jXkoSI3dBgaD\ngcLCQst6UVERBoNhOC8phBAD0thzA7oVhvX/pMTERD788EMA0tPT8fLy6tEMI4QQI03RKlYt9sqm\nGntycjIpKSlUVVURHBzMzp076ejoALrGJF6zZg379+8nNDQUV1dX3n///SEJerQtX758tEO4YxLz\n8LO3eME+Yx4KGp0dt7NYYUyMx64oypC1nQkhxjdb84WiKFxdMtuqsqGpF+0yN9lxF3whhBgcZZy3\nsUtiF0JMOPbcfm4NSexCiAlnvPeKkcQuhBgVbW1t1NfX09DQQENDA01NTTQ1NdHS0kJrayttbW20\nt7fT3t5OZ2cnJpNpyNq7bWmKKSws5IknnqCiogJFUXj66ad54YUXepR74YUXOHDgAHq9ng8++MAy\n8fXBgwfZsWMHJpOJp556ipdeemnQsfRFErsQwmaqqlJfX09JSQllZWWWpaKigsrKSqqqqqiurqam\npoaamhrq6upobW0dtXg1DtpBH+vg4MCvfvUroqOjaWxs5O677yY+Pp7w8HBLmR+Ok3XixAmeeeYZ\n0tPTMZlMPPfccxw6dAiDwcCCBQtITEzsduxQkMQuhBhQe3s7+fn55OXlkZeXR0FBAQUFBRQWFlJU\nVERxcTHNzc13dlJFg85Rj85Bj9bBGa2DMxqdExqdIxqtIxqNDo1Wh6JoUTRaFEVDY3UeN8sv23w/\ntrSxT548mcmTJwPg5uZGeHg4JSUl3ZJzb+NklZWVcePGDUJDQwkJCQEgKSmJL774QhK7EGJ4NDc3\nk5ubS25uLlevXrUs169fp6ioaMBmEK3OGUdXbxz13jjqvXB08cLRxRMHZw8cnN2/WzzQOujRObii\n0TlZN06UqnYtQOHFL4cmsffRFJNWVc/x6garz5OXl0dmZiaxsbHdtvc2TlZxcTElJSU9tp84ceIO\nox+YJHYhJpiqqiouXbpEVlYWV65csSwFBQV9H6RocHb3w8XDH2cPf5zd/XBy9cXJ3Q9nN18cXbzR\nObp2P0btfbAu1azCKHcN7yuxL/H3ZIm/p2X9/2WX9HmOxsZGNm3axFtvvYWbm1uP/aPZ/10SuxDj\nVFNTE5cuXeL8+fNcuHCBixcvcvHiRSoqKnotr2i0uHhORu9jQO8VhN4rEGfPybh4BuDs7o9G0Vpq\nzrcn7bGQrO+Erd0dOzo6ePDBB3n88cdZv359j/29jZNlNBrp6Ojotr2wsBCj0WhTLL2RxC7EOFBa\nWkpmZiZnz54lMzOTc+fOcfXq1V5rjQ46J/R+IbhOMuLmOwVXXyN6HyPOHgEoivb7BH17jdtsR5l7\nALZ0d1RVla1btxIREcGOHTt6LZOYmMiuXbtISkrqNk6Wr68vubm55OXlERQUxJ49e9i9e/egY+mL\nJHYh7IiqqhQXF3P69GnLcubMGcrKynqU1QIGHJmCE0acMOKIl5eB66uTKJkW3r3GrapdTdnjKHn3\nx5bujqmpqXz00UfMnTvX0oXx1VdftTRl9TdOlk6nY9euXaxevRqTycTWrVuH/MUpjKHE3traik6n\nQ2fPEw0KMcSqq6s5efIkGRkZnDx5kpMnT/Y6C5krGu7SOBGidC1TcMKgOqFFsbx7VFWo1zpzQya2\nQeMw+EHAli5ditlsHrDcrl27et2ekJBAQkLCoK9vjTGTRV1cXADQarW4urri4eGBt7c3kyZNwt/f\nn6CgIIKDg5kyZQrTpk0jNDQUd3f3UY5aiKHT3t7O+fPnSU9PJz09nRMnTnD16tUe5dwUDaFaZ0I1\nTkzXOjNN40QADoCCydSVwM3mrhYVOxy/akTIkAIjSsFkMlFfX099fT1FRUX9lg4KCiI8PJzZs2db\nfiyKjIzE0dFxhOIVYvDKyso4fvw4x48fJy0tjdOnT/f4aMcRhZkOToQ5ODNT58wMrTOBigOqqmA2\ndzV7/7BGLqyj2PP0SFYYM4nd3TOC0FnPoppNmMxtmE2tdHY209nZSGdHAx0ddXS019HWVk17WxVt\nrZWUlJRQUlLCP//5T8t5HB0diYqKYtGiRSxevJilS5d26zcqxGgwm81kZWWRmppqWa5fv96j3BQH\nB2Y7uRDh6Ey4gzNTNd83p5jN39XEJYnbTEZ3HGGKRotOowedHkcnnz7Lqai0t9fQ2lJKS3MJLc1F\ntDQV0tZaYWmL/PWvfw3A1KlTWbZsGStWrOBHP/oRU6ZMGanbERNUW1sbJ0+e5NixYxw7dozU1FTq\n6uq6lXHRKMxxdmauiwtzvkvmHhqtpSZuNmNpWhFDSyNNMWOTomhwcp6Ek/MkPL3nfLcRzKZWmhrz\naWq4TlPjDZrqr5Gfn8+HH35omaZv1qxZrFq1ivvvv5/ly5ej1+tH8U7EeFBfX09aWhpHjx7l6NGj\nZGRk0NbW1q3MZAcd81xdiNE7E+XiQqiTE5pbTSo/qI2L4Sc1djuj1bng4RWGh1cYikZBxUxrcwkN\nN3NoqMumoS6H7OxssrOz+fWvf42zszMrVqzgxz/+MevWrZNmG2GV6upqjh07RkpKCkeOHCEzM7NH\nT4mZeifmu+uZ56bnblcXAnU6zGYV1axakvgE6V045sjLUzunKBr0bsHo3YMJMN6Hipmmxjzqqy9R\nV3OR5vo8Dhw4wIEDB3j22WeJiYlh/fr1bNy4kcjISOvGshDjXnl5OUeOHCElJYWUlBQuXrzYbb9O\nUYjy1LPAy40FnnrudtfjqdGimsyYTV3JXDVJFh8rFBu6O9qDcZ/Yb6coWtw8Q3H3moEhdD2dHfXU\nVV+krvIsN6sukZmZSWZmJi+//DIzZszgwQcfZNOmTcybN0+S/ARSXFxsSeIpKSlkZ2d32++s1RDj\n7coiX3cWebsR7anHRdGgmlTMZjOqSRL5WCZNMeOcg5MnfoYl+BmWYDZ3UF9zmdryM9RWnCU3N5fX\nXnuN1157jWnTpvHwww/zyCOPEBUVJUl+nCkoKCAlJYXDhw+TkpLCtWvXuu3X6zQs8PPgHj8PFvu5\nE+XpiqOCJYGbpUZuX7RSY58wNFoHvP2j8A6IAszU1+ZSU3qSmtLTXL9+3ZLkZ82axSOPPEJycjJh\nYWGjHba4Q6qqcuPGDUsiP3LkCHl5ed3KuDloWRzoyZLJntwT4EmUjys6VelqVrE0rwz89aEYo6TG\nPjEpGi2ek8Lx9I/grugnaKjJpro4g+qik2RnZ/PKK6/wyiuvEB0dTXJyMklJSdKNcoxSVZXs7GxL\nG/mRI0d6fPzm6aTjHqM3S4O8WBLkxVwfPVpV6aqRm6WNfLyRl6cCRdHg6R+BV0Ak0+Zt5mblZSoL\n0qguPMnZs2c5e/YsL730EnFxcTz66KNs2rSJSZMmjXbYE5bJZOL8+fMcPXqUI0eOcPTo0R5D1fq4\nOLAkxJe4KT4sDfZmtq8bWhXUDjOq6dYiiXzckhq7+CFFo8Vr8hy8AucQuvBJakvPU5mfSk3hGUsf\n5ueff57Vq1eTnJzMAw880Osg/GLotLS0WD4GOnr0KGlpadTX13crE+DuxNLpk4ibNoklU32I8HVD\nMatg6lpUk7nr92JikBq76ItG68ikKQvwu2shps5WqgpOUnkjldqic+zbt499+/bh4uLCunXrSE5O\n5v7778fZ2Xm0w7Z75eXlpKWlkZqaSlpaGqdOnaKjo6NbmRA/V+Jm+rN0+iSWTJtEqK/rd0n8u5p4\np3nCDFErelJsmMzaHkhiHyI6Bxcmz1hGYNgKOtrrqbyeRnnuMW6WXubPf/4zf/7zn/Hw8OCBBx7g\n4YcfJj4+Hicnp9EOe8zr7OzkwoULlsGyjh8/3qPHiqLA3BAflob5c89MP5bO8MPg5QKdqqVpBZOK\nak9T/IhhJd0dxR1zdPHEOHcNwVFraWmsoOJqGuXZR6ivuM6f/vQn/vSnP+Hp6UliYiIPPvggq1at\nsgxbPJGpqkpeXp5l/PETJ05w+vRpWlpaupVzdXYgNjyAxWGTuSfcn8Uz/PB00YHJhNp5q1YuPVZE\nP6QpRtjCxd2fkPkbuSv2QZpvllKRfZSyK8e4WXHDkuT1ej33338/DzzwAGvXrsXX13e0wx52t7oc\nZmZmcubMGU6fPs2pU6eorq7uUTY02JvYyCAWRQSyOGIyc0J80aFCp6lrlCyTqWugFcvJR/BGBAAK\nin39RCSJXQwVVx8DoUuTmbXscZrriim7coySrGPUFefw17/+lb/+9a9oNBoWLVrEmjVruP/++4mJ\niUFj52NHNzY2kpWVxYULFzh37pxluXnzZo+yk7xdWTjHyPzZBhZEBLEwPJBJHk7Q2Qkm8/dJ3GQa\nhTsR44U0xYhh4eZrZOa9yYQtf5S2hkpKr6RTcjmNiquZpKWlkZaWxs9+9jN8fHxYsWIFK1asIC4u\njsjISLTasffiR1VVKioqyMnJIScnh8uXL3P58mWysrJ6fPxzi/8kd2JmB3P3bCPREQYWzDEQ7OeG\noqrf1cTN39fIhRhKUmMXw03vFcDMpRsIX/YgpvYWyq6eovRKBsVZGdRUl/KXv/yFv/zlLwB4eHiw\ncOFCFi5cSExMDFFRUUybNm3Yk72qqlRXV1NUVERhYSH5+fnk5+dz48YNrl27xrVr12hoaOj1WAcH\nHeGzgpgdZmROhJG54QaiIwwE+rmD2QRmE6qp61c6O2XsWjH8bPgp+Mknn2Tfvn34+/tz4cKFHvvf\neOMNPv74Y6Dr5f/ly5epqqrCy8uLkJAQPDw80Gq1ODg4kJGRMeg4+iOJfYxxcNYzNXoZ0+evQKNR\naKwpoSwnk5IrpynNOUd9ZQmHDh3i0KFDlmOcnJyYMWMG06dPZ+rUqRiNRgICAvD19cXDwwO9Xo+D\ngwOKomA2m+ns7KSjo4O2tjZaWlpoamqisbGR+vp66urqqK2tpbq6mqqqKiorKykrK6OsrKxHl8Lb\neXq6MnOmkZmhBmbNNBA+y0DEzCBC7/LHQYsliXdbhBgFtozuuGXLFp5//nmeeOKJXve/+OKLvPji\niwB8+eWXvPnmm3h5eXVdV1E4fPgwPj59TyI0FCSxj2GKouAZEIxv0FTmrtyITqeh5WYlFdcuUX49\ni/IbVyi7nk19VTkXL17sMZTsUPP08iAoKADjlMkEGwMICZlMSEgA06dPJnR6AJN83FDoJXmbTWDu\nHNbYhLgjNjTFxMXF9dm8eLtPPvmE5OTkbtvUEfiJVBK7nXH3DcDbP5DIuHh0Wg0OWg2dbc3UlRVS\nW1pAfVUpN6sqaKyrprmulraWJtpbWzB1dqIAGo0GnU6HTqfDyckJJ2dn9HpXXN3ccPdwx8vLE09P\nD3x8vfD19cI/wAc/fy+CJvui1+vQKJ1oFBNaOr/7fSeK0olWMaGoJumRIuxCXy9PU66Uk5Jd0eu+\nO9Xc3Mw//vEPfvvb335/XUVh5cqVaLVatm/fzrZt24bkWreTxD4OOLu6EzxrNtMi5+KgU9DptDjo\nNDjpNDhqNeg0Cg5aDQ5aBQdN17pOUdBqupbvfw86RUWjmC0J/Fby1tIJSK1bjBN9DNu7LDKQZZGB\nlvX/2jv4n4L//ve/s3TpUkszDEBqaiqBgYFUVlYSHx9PWFgYcXFxg75GX2zuR3fw4EHCwsKYMWMG\nr7/+eo/9hw8fxtPTk5iYGGJiYvjlL39p6yWFEMI2WsW6xQaffvppj2aYwMCu/zT8/PzYsGHD2Hx5\najKZeO655zh06BAGg4EFCxaQmJhIeHh4t3LLli1j7969NgUqhBBDZbj7sd+8eZMjR47wySefWLY1\nNzdjMplwd3enqamJr776ipdffnlYrm9TYs/IyCA0NJSQkBAAkpKS+OKLL3ok9pF4WSCEEFazoVdM\ncnIyKSkpVFVVERwczM6dOy09xrZv3w7A3/72N1avXt1tqJDy8nI2bNgAdHWDfOyxx1i1apUNN9E3\nmxJ7cXExwcHBlnWj0ciJEye6lVEUhbS0NKKiojAYDLzxxhtERET0OFdbayWlxfsAcHOfgbvHTFtC\nE0KME3Wll6grvQRAfUXOkJzTlok2du/ePWCZzZs3s3nz5m7b7rrrLs6ePTvo694JmxK7NfN+zps3\nj8LCQvR6PQcOHGD9+vXk5PT8w3Fy9iPQsNaWcIQQ45BXYCRek7sqg4UXv6S2+JztJ7XzYToGYtPd\nGQwGCgsLLeuFhYUYjcZuZdzd3dHr9QAkJCTQ0dFBTU2NLZcVQgjbaDTWLXbKpsjnz59Pbm4ueXl5\ntLe3s2fPHhITE7uVKS8vt7SxZ2RkoKrqsH91JYQQ/Rrnid2mphidTseuXbtYvXo1JpOJrVu3Eh4e\nzrvvvgt0vUj47LPPeOedd9DpdOj1ej799NMhCVwIIQZNRnfsX0JCAgkJCd223XozDPDss8/y7LPP\n2noZIYQYOnZcG7eGfHkqhJh4dOM79Y3vuxNCiN5IjV0IIcYZSexCCDHOSGIXQohxRnrFCNE1zPrg\n/ynIWEHjnWpvf8bjvMY+vu9ODBnb6jfju3YkQLG3P2M7+EDJZDJZpti7U1JjF1aRGrvoj73V2BXt\n2E99Wq2WY8eOoaqqVeNy/dDYvzshhBhqdtIUEx0dzQMPPMBDDz1kGXNLURQ2btzY73GS2IUQE4+d\nvDxtbW3Fx8eHb775ptt2SexCCHE7O6mxf/DBB4M6zj7uTgghhpIdvDwFyM7O5r777iMyMhKA8+fP\nWzVv9OhHLoQQI81OEvu2bdt49dVXcXR0BGDOnDlWzeAkTTFCiIlnDCRtazQ3NxMbG2tZVxQFBweH\nAY+zj7sTQoihpNNat/TiySefJCAggDlz5vS6//Dhw3h6ehITE0NMTEy3ppODBw8SFhbGjBkzeP31\n1wcM08/Pj6tXr1rWP/vsMwIDAwe+vQFLCCHEeGNDjX3Lli08//zzPPHEE32WWbZsGXv37u22zWQy\n8dxzz3Ho0CEMBgMLFiwgMTGR8PDwPs+za9cunn76aa5cuUJQUBB33XUXH3/88YAxSmIXQkw8NiT2\nuLg48vLy+i1zazrQH8rIyCA0NJSQkBAAkpKS+OKLL/pN7NOnT+ef//wnTU1NmM1m3N3drYpRErsQ\nYuLpI7EfTr9KSvrVXvdZS1EU0tLSiIqKwmAw8MYbbxAREUFxcTHBwcGWckajkRMnTvR7rqqqKnbu\n3MmxY8dQFIW4uDh+8Ytf4Ovr2+9xktiFEBOP0ntiX754JssXz7Ss73zz4B2fet68eRQWFqLX6zlw\n4ADr168nJydnUGEmJSWxbNky/vrXv6KqKp988gmPPPIIhw4d6vc4eXkqhJh4FI11yyC4u7tbPv9P\nSEigo6ODmpoajEYjhYWFlnKFhYUYjcZ+z1VWVsbPf/5z7rrrLqZNm8bPfvYzysvLB4xBErsQYuIZ\nxn7s5eXlljb2jIwMVFXFx8eH+fPnk5ubS15eHu3t7ezZs4fExMR+z7Vq1Sp2796N2WzGbDazZ88e\nVq1aNWAM0hQjhJh4NL13ZbRGcnIyKSkpVFVVERwczM6dO+no6ABg+/btfPbZZ7zzzjvodDr0ej2f\nfvopADqdjl27drF69WpMJhNbt27t88Wpm5ubZUTHN998k5/85CcAmM1mXF1d+d///d9+Y5TELoSY\neAbZzAIM+OXns88+y7PPPtvrvoSEBBISEga8RmNj46Biu0USuxBi4rGTL0+ha3yYvLw8Ojs7Ldtk\ndEchhLidDTX2kbRlyxYuXLhAZGQkmh/8ZySJXQghbmcnif3EiRNcunTpjmdQso+7E0KIoaQo1i2j\nbMGCBWRlZd3xcVJjF1aRyaxFfxQU+5r31A7mPIWuppjFixczefJknJycgK4vW8+fP9/vcfZxd0II\nMYQUO2mK2bp1Kx999BGzZ8/u1sY+EEnsQoiJx04Su7+//4AfMfVGErsQYuKxk8QeExPDo48+yrp1\n6yyzKCmKIr1ihBCihzHwYtQazc3NODo68tVXX3XbLoldCCFuZycfKH3wwQeDOk4SuxBi4rGTppgt\nW7Z0W7/Vn/0Pf/hDv8dJYhdCTDwa+0h9a9eutSTzlpYWPv/8c4KCggY8zua7O3jwIDt27MBkMvHU\nU0/x0ksv9SjzwgsvcODAAfR6PR988AExMTG2XlYIIQbPTppiNm3a1G390UcfZcmSJQMeZ9Pd3Zqc\n9eDBg2RlZbF7924uX77crcz+/fu5evUqubm5vPfeezzzzDO2XFIIIWw3jBNtDKecnBwqKysHLGdT\njd2ayVn37t3L5s2bAYiNjaWuro7y8nICAgJsubQQQgyenfSK+eG47IqiEBAQwOuvvz7gcTYldmsm\nZ+2tTFFRUY/E3tZaSWnxPgDc3Gfg7jETIYSoK71EXeklAOorBjd3aA9jsDbem8GOy25TYrd2xLFb\n00T1d5yTsx+BhrW2hCOEGIe8AiPxmhwBQOHFL6ktPmf7Se0ksUNX5Tg/P7/beOz33ntvv8fYlNgN\nBsOAk7PeXqaoqAiDwWDLZYUQwiaqnST2l156iT179hAREYFW+/10fsOa2H84OWtQUBB79uzpMW1U\nYmIiu3btIikpifT0dLy8vKR9XQgxqsyqw2iHYJXPP/+c7Oxsy8iO1rIpsfc1Oeu7774LdE3sumbN\nGvbv309oaCiurq68//77tlxSCCFspqqDr7E/+eST7Nu3D39/fy5cuNBj/8cff8z//M//oKoq7u7u\nvPPOO8ydOxeAkJAQPDw80Gq1ODg4kJGR0e+1pk+fTnt7+8gmduh9ctbt27d3W9+1a5etlxFCiCGj\n2tDTe8uWLTz//PM88cQTve6fNm0aR44cwdPTk4MHD/L000+Tnp4OdL1fPHz4MD4+PlZdy8XFhejo\naO67775u47G//fbb/R5nH59fCSHEEDLbUGOPi4sjLy+vz/2LFy+2/D42NpaioqJu+2/vTNKfxMRE\nEhMTLR1OVFW1qtOKJHYhxISjou11+5GUUxxJOT1k1/n973/PmjVrLOuKorBy5Uq0Wi3bt29n27Zt\n/R7/L//yL4O6riR2IcSE01cbe9y9C4m7d6Fl/b9feXfQ1/j222/5wx/+QGpqqmVbamoqgYGBVFZW\nEh8fT1hYGHFxcYO+Rl/so8+PEEIMIRWNVctgnT9/nm3btrF37168vb0t2wMDAwHw8/Njw4YNA748\nHSxJ7EKICcek6qxaBqOgoICNGzfy0UcfERoaatne3NxMQ0MDAE1NTXz11VfMmTNnSO7ndtIUI4SY\ncPpqY7dGcnIyKSkpVFVVERwczM6dO+no6AC6egS+8sor1NbWWgY8vNWtsayszDLzUWdnJ4899hir\nVq3q91rZ2dm88cYb5OXlWb48VRSFb775pt/jJLEPt7Ew1tCgYuh+kDrY01iOHkFj4ZlPMOpI/xnb\nyJZ+7Ld/hHm73/3ud/zud7/rsX3atGmcPXv2jq710EMP8cwzz/DUU09ZvjyVXjFCCNELW9rPR5KD\ng8OghjqXxC6EmHBs6cc+EmpqalBVlXXr1vGb3/yGjRs3dvv6dKAPnCSxCyEmHFva2EfCvHnzujW5\nvPHGG93237hxo9/jJbELISYcW9rYR0J/X7ZaY2zfnRBCDAOzqrNqGW2/+c1vqK2ttazX1tby29/+\ndsDjJLELISYcM1qrltH23nvvdfvAydvbm/fee2/A40b/vyQhhBhhY70p5haz2YzZbEaj6YrXZDJZ\n+sz3RxLjmlkkAAAa3klEQVS7EGLCsZfujqtXryYpKYnt27ejqirvvvsu999//4DHSWIXQkw49lJj\nf/3113nvvfd45513AIiPj+epp54a8DhJ7EKICWcstJ9bQ6vV8swzz9zxR0qS2IUQE85Y6PFijZyc\nHP7jP/6DrKwsWlpagK4hBa5fv97vcfbx84gQQgwhs6patYy2LVu28K//+q/odDq+/fZbNm/ezGOP\nPTbgcZLYhRATjlm1bhltLS0trFy5ElVVCQkJ4T//8z/Zt2/fgMfZx88jQggxhMZCbdwazs7OmEwm\nQkND2bVrF0FBQTQ1NQ14nCR2IcSEMxZq49Z48803aW5u5u233+bnP/859fX1/PGPfxzwOEnsQogJ\nx15q7AsXds2/6u7uzgcffGD1cZLYhRATzlhP7OvWrUNRFNRe4lQUhb179/Z7vCR2IcSE0zHG22LS\n09MxGo0kJycTGxsLYEnyMoOSEEL0YozndUpLS/n666/ZvXs3u3fvZu3atSQnJxMZGWnV8dLdUQgx\n4djSj/3JJ58kICCAOXPm9Hn+F154gRkzZhAVFUVmZqZl+8GDBwkLC2PGjBm8/vrrfR6v0+lISEjg\nww8/JD09ndDQUJYtW8auXbusuj9J7EKICceWfuxbtmzh4MGDfZ57//79XL16ldzcXN577z3LcAAm\nk4nnnnuOgwcPkpWVxe7du7l8+XKf52ltbeUvf/kLjz/+OL/5zW/493//dzZs2GDV/UlTzHAbCz/y\nDSqG7gcN3KrXH9uOvmNj4ZlPMAoKqh09+L5q42fSjnIm7Vi/x8bFxfU7w9HevXvZvHkzALGxsdTV\n1VFWVsaNGzcIDQ0lJCQEgKSkJL744gvCw8N7nOMnP/kJly5dYs2aNfziF7/o96eD3khiF0JMOH0l\n9ujFS4levNSy/vv/fe2Oz11cXExwcLBl3Wg0UlxcTElJSY/tJ06c6PUcH3/8Ma6urrz11lu89dZb\n3fYpikJ9fX2/MUhiF0JMOMP98rS3bop3wmw223S8JHYhxITTOYyZ3WAwUFhYaFkvKirCaDTS0dHR\nbXthYSFGo3FYYpCXp0KICWc4R3dMTEzkww8/BLr6o3t5eREQEMD8+fPJzc0lLy+P9vZ29uzZQ2Ji\n4lDeloXU2IUQE44tX54mJyeTkpJCVVUVwcHB7Ny50zIP6fbt21mzZg379+8nNDQUV1dX3n//faCr\nC+OuXbtYvXo1JpOJrVu39vridChIYhdCTDi2tMTs3r17wDJ99TdPSEggISFh8Be3kiR2IcSEYxrj\nY8XYatCJvaamhkceeYT8/HxCQkL485//jJeXV49yISEheHh4oNVqcXBwICMjw6aAhRDCVmN9SAFb\nDfrl6WuvvUZ8fDw5OTncd999vPZa7/09FUXh8OHDZGZmSlIXQowJ9jI13mANOrH/8OuqzZs387e/\n/a3Psrb26RRCiKHUaVatWuzVoJtiysvLCQgIACAgIIDy8vJeyymKwsqVK9FqtWzfvp1t27b1Wq6t\ntZLS4q65/NzcZ+DuMXOwoQkhxpG60kvUlV4CoL4iZ0jOac+1cWv0m9jj4+MpKyvrsf2///u/u60r\nitLnGMGpqakEBgZSWVlJfHw8YWFhxMXF9Sjn5OxHoGHtncQuhJgAvAIj8ZocAUDhxS+pLT5n8zlt\n/LBzzOs3sX/99dd97gsICKCsrIzJkydTWlqKv79/r+UCAwMB8PPzY8OGDWRkZPSa2IUQYqSY7biZ\nxRqDbmNPTEy0TKr6xz/+kfXr1/co09zcTENDAwBNTU189dVXdzxKmRBCDDV5edqHn/70p3z99dfM\nnDmTb775hp/+9KcAlJSUsHZtV5NKWVkZcXFxREdHExsby49//GNWrVo1NJELIcQgmcyqVYu9GvTL\nUx8fHw4dOtRje1BQEPv2db0EnTZtGmfPnh18dOKOtTY1Ul9dQUt9La0NdXS2NNPR1oJqMqFRVBy+\n+57AycEBJycnXPQuuOpd8XB3x8PTHS9PD7y83PH01I/4MOpCjJTx3hQjX57aIVVVqSsvoTIvm8r8\nq1QWXKOqOI/qkgJaGvofp9laOp0Obx8fJvn74u8/iYCASUwO8sNg8Mdo9GPqFH+mTvXDb5IbVsyt\nK8SY0mka329PJbHbgc72NkpzsijLPUdp7jlKcy7Q0lDXa1kXFxcMBgMBAQH4+vri4eGBXq/HwcEB\nRVEwm810dnbS0dFBW1sbLS0tNDU10djYSH19PTdv3qS2tpbGxkYqKyqorKig78m7wM1Nz13Tggid\nHsSM0CBCQwOZNcPArNBA/HzdpNIvxiSpsYsRp5rNVOXnUHI5g+Irpyi/eh5TR3u3Mr6+vsybN4+o\nqCgiIyOZNWsW06dPx8/Pr8+up3eira2N6upqKioqKCsro7S0lNLSUgoLCykqKiI/P5+8vDwaGhq4\ncP4qF85f7XEOHx93IsKMRH63zA4LYm64AR9PF5vjE8IWktjFiOhobaIo9xQlWccpuZxOa2P3Gvnc\nuXO59957WbJkCYsWLWLq1KlDksD74uTkRFBQEEFBQX2WUVWVmpoabty4YZm8Nzc3l+zsbK5cuUJN\nTT3H0i5zLK17nd8Y5EN0pLFriTAwL9LAlEBPqd2LEWPPPV6sIYl9FLU21lKefZySrGNUXsvEbOq0\n7JsyZQqrV69m1apVrFixAl9f31GMtHeKouDr64uvry/z58/vtk9VVUpLS7l06RKXLl3iwoULXLx4\nkYsXL1JUUkNRSQ1ffn3eUt7HS8/8OcHcHRnE/EgDCyODCPJ1G+lbEhOEPfd4sYYk9hHW1lhDRc5x\nSq8cpTr/AnxXc9BoNCxdupR169bx4x//mPDw8GGtkQ83RVEsNf74+HjLdpPJxLVr1zh79iyZmZlk\nZmZy+vRpqqqq+OpoNl8dzbaUNfh7EDvbQOzsIGLDJ3P3DH/0DjLp11ikYl+JUppihM3am29Sce04\n5TlHqS28aEnmjo6OxMfHs2HDBhITE/Hz8xvlSIefVqtl5syZzJw5k4cffhjoqt0XFhZy6tQpTp06\nRUZGBidPnqS4op6/flPPX7/pasrRaTVEhfqxODKQe8Inc0+YP8E++tG8HWGnpFeMGJTO9maqrp+g\nPPcYNQVnUdWuv0iOjo7cf//9PPTQQ6xbtw5PT89RjnT0KYrClClTmDJlChs3bgS6ZmnPycnhxIkT\nHD9+nPT0dC5cuMDp7HJOZ5dza34a4yRXloRPZmm4P0tm+REZ5Il29G5F2AmpsQurmTvbqSk4R8W1\no1Tnn8Js6poHUavVsmpVAklJSTzwwAOSzK2g0WgICwsjLCzMMjx0Y2MjGRkZHD9+nNTUVI4fP05R\nVR17jl5jz9FrAHjqHVg8w4+lM/1YGurH3UZvHDX226Qlhoetif3gwYPs2LEDk8nEU089xUsvvdRt\n/xtvvMHHH38MQGdnJ5cvX6aqqgovL68RmXxIEruNVLOZuvIsKm+kUpV/gs72Zsu+e++9l+TkZDZt\n2sSkSZNGMcrxwc3NjR/96Ef86Ec/Arpq9VlZWaSmpnL06FGOHTtGfn4+B8+VcPBcCQDODloWhviw\ndPok4qZNYqHRC71O6vQTnS1T45lMJp577jkOHTqEwWBgwYIFJCYmdpuY+sUXX+TFF18E4Msvv+TN\nN9+0zDB3a/IhHx8f226iH5LYB0FVVRqqr1JdcJzK/ON0tN607IuJiSE5OZmkpCSCg4NHMcrxT6PR\nMHv2bGbPns327dsBKCoq4siRIxw9epSjR49y6dIljuRWciS3EgCdRmGe0YulIb4sneLDPUYvPHTy\nz2CisWXY3oyMDEJDQwkJCQEgKSmJL774olti/6FPPvmE5OTkbtuGe/Ih+RttJVVVaarLp6roBFWF\nx2lrqrLsCw0NJTk5mUcffZSwsLBRjFIYjUYeffRRHn30UQCqq6stST4lJaVrisaCWjIKavl/gEaB\n6MmeLA32ZqnRm8UBnng5SI1+vOurKebGuQxunO+/aaS4uLhbpc1oNHLixIleyzY3N/OPf/yD3/72\nt5Zt1k4+ZAtJ7P1QVZXmm4XUlJ6kujiD1sbvJx0JCgoiKSmJ5ORk7r77brvumjie+fr6sn79esuw\n0jdv3iQtLY2UlBRSUlI4deoUZ0pvcqb0Jm9n5KEAc/zcWBLoxdJALxb7e+DtKP9Mxpu+EvvUOQuY\nOmeBZf3bj3b1KHMn/9b//ve/s3TpUkszDFg/+ZAt5G/sbVTVTNPNPGrLz1BTdprWpu+n/PPz82PT\npk0kJSWxdOlSNBrpU21vPD09SUhIICEhAeiaJ+D48eOkpKRw+PBhMjIyOF/ZyPnKRt45XwRAhLee\newI8ucffk9hJ7vhJord7tnR3NBgMFBYWWtYLCwsxGo29lv300097NMOMxORD8jcUMJs7qK+5Ql3F\nOWorMulo+/5zfn9/fzZs2MDDDz/Mvffei07aY8cVV1dXVq5cycqVKwFoaWkhPT3dkujT09PJqm0m\nq7aZ310pBSDU3ZnFfh4s8nUn1seNQEeH0bwFMQi2fHk6f/58cnNzycvLIygoiD179rB79+4e5W7e\nvMmRI0f45JNPLNuam5sxmUy4u7tbJh96+eWXBx1LXyZslmpvreFmzUVuVl/gZnUWZlObZZ/RaGTj\nxo1s2LCBuLg4tFppc50oXFxcWLFiBStWrAC6BkPLyMiwJPrjx49ztaGZqw2t/Ol6BQBT9I7E+riz\n0MuVhZ5uBDtN2H9WdsOW7o46nY5du3axevVqTCYTW7duJTw8nHfffRfA8iL/b3/7G6tXr8bF5ftB\n78rLy9mwYQPQ1Q3yscceG5bJhxR1uF/PWhOEouDuGUHorGfv5KCeE0Eo3du/FI1i2WYytdBwM5f6\nuivU116mtamk26FRUVGsW7eO9evXM2/ePGkzF71qb2/nzJkzpKSkcOTIEY4dO0Z9ffcx8P0cdSzw\ndGW+u5673VyY6eyERgXVpGI2q6hmFbO56wNks7n7Ytn23a8mU9e2H5btdpz6/X5V/b78rXK3L/W+\nAZyK30DxXWHdemaoqtpVxqxiGR1Ava33xq2L9bKv23H0cXx/Zfujfn/dwotfcj3jQ5t6lSiKwv/3\n+fmBCwL/s2HusPdgGQ7jtmrR0dFAc+N1Guuv0nAzl+bGAn74t8nNzY377ruP+++/n7Vr10rXRGEV\nR0dHFi1axKJFi3jppZcwmUycP3+elJQUjh49ypEjR6isqmJ/5U32V3Z1g3XVaohxdWGeqwvRemfm\nuDjjPPhZKcUQkC9P7UBnRyMtzcU0NxfS0phPU2M+7W1V3cpotVoWLlzIfffdR3x8PIsWLcLR0XGU\nIhbjhVarJSYmhpiYGHbs2IGqqmRnZ1u6WB47dowbN25wrL6JY/VNXccAs5ydiNG7MNfFhdlOzvhr\ndMhchCNHEvsYoKoqJlMLHe21tLfV0NZWSVtrJW2t5bS0lNLZ0XM6OFdXVxYsWEBcXBxLly7lnnvu\nwc1NhoEVw0tRFMtQCLf6J5eUlJCamsqxY8c4duwY586dI6u1jazWNj6m60W9v07HbCdnZjs5E+no\nwgwHJxykVj9sJLGPkNaWUgpu7EZVOzCbOzCZWjB1NtPZ2UhHRz2quaPPY11dXZkzZw7R0dHMnz+f\n+fPnExkZKT1YxJgQFBTEQw89xEMPPQR8P+ZNWloaqamppKenU1FXxzedjXzT1Ah01epDHZwId3Am\nzMGZWQ7OGBRHpFY/NEwdptEOYViNmczX0V5LdeWxPve7uroSHBzM1KlTmTZtGqGhocyaNYvw8HBC\nQkKkT7mwG72NeZOdnU16erplJMtLly6R3dFGdkcb8F1bvaIhVOfETK0zoVpnpitO+CsOSLK/c1Jj\nHyFvv/02Wq0WFxcX9Ho9Hh4eeHt7M2nSJAICAnB3dx/tEIUYFhqNhvDwcMLDw9myZQsADQ0NlrHp\nT5w4wcmTJykqKuJcRwvnOlosx7qiYbrWiWkaZ+7SODFVdSIQRxRJ9v2SxD5Cnn/++dEOQYgxw93d\nvVt/eoDS0lJOnTrFyZMnOX36NKdPn6a8vJzzphbOm75P9o4oTFGcmIIjU3HCiBPBqhOuMlK9hSR2\nIcSYEBgYyLp161i3bh3Q1amgpKSEM2fOcObMGct0g/n5+VxVW7lKa7fjvdExub4CbaYG6vNx9Q3G\n1TcYB+eJ16nAbJLELoQYgxRFwWAwYDAYLMkeoLa2lnPnznHu3DkuXLjQ1QsnK4va5mZqO+ogO6Vr\n+Y6jqzeuPkb0Pgb03gb0XkHovYNwcQ9A0Y7PFCE1diGEXfH29mb58uUsX77css1sNnPjxg0uXbrE\nxYsXycrKIisriytXrtDSVEt7Uy21hRduO5OCs4cfLp6TcfEMwMUjAGd3f5w9/HF288PRxRNFY5/N\nO6okdiGEvdNoNEyfPp3p06eTmJho2W42myksLOTKlSvk5OSQnZ1Nbm4uOTk5FBQU0FpfQWt9BbWF\nPc+paLQ4uvrg7OqLk5svjnofnFy9cdR/tzh74uDsgYOz+5j7D6CzU7o7CiHGKY1Gw9SpU5k6dSqr\nV6/utq+9vZ38/HyuXbvGjRs3uH79Ovn5+eTl5VFQUEB5eTltDZW0NVQOcBUFnZNrV5J3ckPn5IrO\n0Q2dox6tox6tgws6B2e0Ds5odE5otI5odY4oGgc0Wh2KokVRNCiKQkdL3QDXso60sQshJiRHR0dm\nzJjBjBkzet3f2tpKcXExRUVFFBcXU1JSQnFxMWVlZZSWllJRUUFFRQU1NTV0tjXS2dZIS69nGnnS\nxi6EEL1wdna2NO/0p7Ozk9raWqqqqqitraWmpoa6ujrq6uqor6+nvr6ehoYGmpubaWpqoqWlhZaW\nFtra2mhra6OjowOTyYTJZMJsNnPlyhWbY5fELoQQNtDpdPj5+eHn5zck5xuKIbXl5akQQowzUmMX\nQohxZry/PJWRswbh8OHDox3CHZOYh5+9xQv2GfNQ6Ow0WbXYK0nsg2CP/xgk5uFnb/GCfcY8FMxm\n1arFXg06sf/f//0fkZGRaLVazpw502e5gwcPEhYWxowZM3j99dcHezkhhBgytib2gfLa4cOH8fT0\ntMyu9ctf/tLqY4fCoNvY58yZw+eff26Zkbs3JpOJ5557jkOHDmEwGFiwYAGJiYmEh4cP9rJCCGEz\n1YY2dmvz2rJly9i7d++gjrXVoBN7WFjYgGUyMjIIDQ0lJCQEgKSkJL744oteb2IoujCNpJ07d452\nCHdMYh5+9hYv2GfMtrKlmcXavKaqPa9xJznRFsPaK6a4uJjg4GDLutFo5MSJEz3K9fYAhBBiuPSV\n2GsKzvcyGFp31uQ1RVFIS0sjKioKg8HAG2+8QUREhNU50Vb9Jvb4+HjKysp6bH/11Ve7DRPaF3ur\nhQshJgbVbO51u7dxNt7G2Zb1a6mf9ChjTV6bN28ehYWF6PV6Dhw4wPr168nJyRl8wHeo38T+9ddf\n23Ryg8FAYeH3w8IVFhZiNBptOqcQQtjK3NF7YreGNXnth1N5JiQk8G//9m/U1NRgNBpHJCcOSXfH\nvppS5s+fT25uLnl5ebS3t7Nnz55uQ4YKIcRoUE2qVUtvrMlr5eXllryYkZGBqqr4+PiMWE4cdGL/\n/PPPCQ4OJj09nbVr15KQkABASUkJa9euBbrGiNi1axerV68mIiKCRx55RHrECCFGnWpWrVp601de\ne/fdd3n33XcB+Oyzz5gzZw7R0dHs2LGDTz/9tN9jh/4GR9CBAwfUWbNmqaGhoeprr73Wa5nnn39e\nDQ0NVefOnaueOXNmJMPr1UAxf/vtt6qHh4caHR2tRkdHq//1X/81ClF22bJli+rv76/Onj27zzJj\n7fkOFPNYer63FBQUqMuXL1cjIiLUyMhI9a233uq13Fh51tbEO9aec0tLi7pw4UI1KipKDQ8PV3/6\n05/2Wm4wzxhQV2z/i1XLCKfIITNiUXd2dqrTp09Xb9y4oba3t6tRUVFqVlZWtzL79u1TExISVFVV\n1fT0dDU2NnakwuuVNTF/++236rp160Ypwu6OHDminjlzps8kOdaer6oOHPNYer63lJaWqpmZmaqq\nqmpDQ4M6c+bMMf132Zp4x+JzbmpqUlVVVTs6OtTY2Fj16NGj3fYP9hkD6vJtn1m12GtiH7EhBX7Y\nf9PBwcHSf/OH9u7dy+bNmwGIjY2lrq6O8vLykQqxB2tihrHTXTMuLg5vb+8+94+15wsDxwxj5/ne\nMnnyZKKjowFwc3MjPDyckpKSbmXG0rO2Jl4Ye89Zr9cDXTM5mUwmfHx8uu235Rnb0sZuD0YssffW\nf7O4uHjAMkVFRSMVYg/WxPzD/qpr1qwhKytrpMO02lh7vtYY6883Ly+PzMxMYmNju20fq8+6r3jH\n4nM2m81ER0cTEBDAihUriIiI6LbflmesdpqsWuzViA3ba22f9ttrDaPZF94e+qveqbH0fK0xlp9v\nY2MjmzZt4q233sLNza3H/rH2rPuLdyw+Z41Gw9mzZ7l58yarV6/m8OHDLF++vFuZwT5je66NW2PE\nauzW9P28vUxRUREGg2GkQuzB2v6qt35kTEhIoKOjg5qamhGN01pj7flaY6w+346ODh588EEef/xx\n1q9f32P/WHvWA8U7Vp8zgKenJ2vXruXUqVPdttvyjG3pFWMPRiyxW9N/MzExkQ8//BCA9PR0vLy8\nCAgIGKkQe7Clv+pYNNaerzXG4vNVVZWtW7cSERHBjh07ei0zlp61NfGOtedcVVVFXV0dAC0tLXz9\n9dfExMR0K2PLMx7vbewj1hTzw/6bJpOJrVu3Wvp+Amzfvp01a9awf/9+QkNDcXV15f333x+p8AYd\n82effcY777yDTqdDr9db+quOhuTkZFJSUqiqqiI4OJidO3fS0dFhiXWsPV8YOOax9HxvSU1N5aOP\nPmLu3LmWZPPqq69SUFAAjL1nbU28Y+05l5aWsnnzZsxmM2azmZ/85Cfcd999Q5Yv7Lk2bg1FHWuv\nwoUQYhgpisI9Gz60qmza50+Mud5C1pA5T4UQE854r7FLYhdCTDi2DAJmDySxCyEmHHt+MWoNSexC\niImnj/HYxwtJ7EKICUdq7EIIMc7Iy1MhhBhnpMYuhBDjzThvY5cPlIQQE8qdDMbm7e09ZsbMuRNS\nYxdCTCgToS47YoOACSGEGBmS2IUQYpyRxC6EEOOMJHYhhBhnJLELIcQ4I4ldCCHGmf8fcW3SaCpt\nGikAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the Mach number-area ratio relation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M_range = np.linspace(0, 4, 201)\n",
      "plt.plot(M_range, fl.A_Astar(M_range))\n",
      "plt.plot(x, A_e / A_c * np.ones_like(x))\n",
      "plt.xlim(0, 3)\n",
      "plt.ylim(0, 4)\n",
      "plt.title(\"Mach number and area ratio\")\n",
      "plt.xlabel(\"M\")\n",
      "plt.ylabel(\"A / A*\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/home/juanlu/.local/lib/python3.3/site-packages/skaero/gasdynamics/isentropic.py:150: RuntimeWarning: divide by zero encountered in true_divide\n",
        "  ((self.gamma + 1) / (2 * (self.gamma - 1))) / M\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "<matplotlib.text.Text at 0x7fe7c3cd85d0>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEVCAYAAAAPRfkLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVPXiBvB3WFyGHRTUAUTDBVxYSlHcwBVcyKWuSymW\neVEzteVW3sokzeqmrZZpi2YWWlZKCpbb4IqkQporLqwiioCAgMDw/f1hzs+RRbbDYWbez/PwxMw5\nnPMejp2Xs45CCCFAREQEwETuAERE1HSwFIiISIulQEREWiwFIiLSYikQEZEWS4GIiLRYCtQgpk+f\njjfeeEPuGBUkJSXBxMQE5eXlckepETc3N+zevVvuGJL4/vvvMWLECLlj0AOwFIyEm5sbmjdvjhs3\nbui87+PjAxMTE6SkpNRr+gqFAgqFol7TIMP5PVZWxk888QR+//13GVNRTbAUjIRCoUDHjh0RERGh\nfe/kyZMoKioyiI1QYygrK5M7Qo1JnbWme168N1b/sBSMyJNPPon169drX3/77beYNm2azv+427dv\nh4+PD2xsbODq6orw8HCdaRw4cAD+/v6ws7ODq6urzvSys7MxevRoWFtbo0+fPrh06VKlOe7+Fbl+\n/Xq0b98erVu3xrJly7TD7z8UpVar4eLion3t5uaG5cuXo2fPnrCyssKMGTOQmZmJ4OBg2NjYYNiw\nYcjNzdWZ59dffw2VSoV27dphxYoV2veFEHj33Xfh7u6OVq1aYeLEicjJydHJ+c0336B9+/YYOnRo\nhWXJzc3F6NGj4ejoCHt7e4wZMwbp6ena4QEBAVi0aBH69+8Pa2trjBgxQmdv7bvvvkP79u3RqlUr\nnd9BZapbN1Vl/eabb+Dp6Ql7e3sEBQXp7BHOnz8frq6usLGxwSOPPIIDBw5UOe/p06dj9uzZGDly\nJCwtLaFWq6vNM3DgQACAra0trK2tERsbi3Xr1mHAgAHacQ4dOoRevXrB1tYWvXv3xuHDh6tdfmok\ngoyCm5ub2LVrl+jSpYs4c+aMKCsrE87OziI5OVkoFAqRnJwshBBCrVaLv//+WwghxIkTJ4STk5PY\nsmWLEEKIpKQkYWVlJTZu3CjKysrEjRs3REJCghBCiNDQUOHg4CD+/PNPUVZWJp544gkxadKkSrNc\nvnxZKBQK8e9//1sUFxeLv/76SzRv3lycPXtWCCHE9OnTxRtvvKEdf+/evcLZ2VlnWfr27SuuXbsm\n0tPThaOjo/Dx8REJCQmiuLhYDB48WISHh+vMa8qUKaKwsFCcPHlStG7dWuzatUsIIcRHH30k+vbt\nK9LT00VJSYkICwsTkydP1vnZ0NBQUVhYKIqLiyssy40bN8Qvv/wiioqKRH5+vnj88cfF2LFjtcMH\nDRok3N3dRWJioigqKhIBAQHi1VdfFUIIcerUKWFpaSn2798vbt++LV544QVhZmYmdu/eXenvrbp1\nc3/WoqIisWXLFuHu7i7Onj0rNBqNWLp0qfD399dOb8OGDSI7O1toNBqxYsUK0aZNG3H79u1K5x0a\nGipsbGzEoUOHhBBCFBcXP/DfikKhEBqNRjuNtWvXiv79+2t/b7a2tmLDhg1Co9GIiIgIYWdnJ27c\nuFHp/KnxsBSMxN1SWLp0qVi4cKGIjo4Ww4cPF2VlZTqlcL/58+eL559/XgghxLJly8T48eMrHW/6\n9Oli5syZ2tdRUVGia9eulY57dwOWnp6ufa93795i06ZN2mm9/vrr2mGVlcIPP/ygfT1hwgQxZ84c\n7etPP/1Uu2G+O69z585ph7/88stixowZQgghunbtqrMRvnLlijA3NxcajUb7s5cvX650OSoTHx8v\n7OzstK8DAgLE22+/rX39+eefi6CgICGEEOHh4doCEkKIW7duiWbNmlVZCve7d91UljUoKEh8/fXX\n2tcajUYolUqRkpJS6fTs7OzEiRMnKh02ffp0ERoaWus8VZXC+vXrhZ+fn87P9+3bV6xbt67aeZD0\nePjIiCgUCkydOhXff/99pYeOAODIkSMIDAyEo6MjbG1tsXr1au3hjtTUVHTs2LHK6Ts5OWm/b9my\nJQoKCqrN06ZNG+33SqXygeNXN697X7do0aLCtO49/OTq6oorV64AAJKTkzFu3DjY2dnBzs4Onp6e\nMDMzQ2ZmZqU/e7/CwkKEhYXBzc0NNjY2GDRoEG7evKnze713Oe/9vVy5cgXOzs46vwMHB4cq51Xd\nuqksa3JyMubPn69dtrvTvnt4a/ny5fD09IStrS3s7Oxw8+ZNZGVlVTn/+38PNclTlStXrsDV1VXn\nvfbt2+sceiN5sBSMjKurKzp27Ijo6GiMHz++wvApU6Zg7NixSEtLQ25uLmbNmqXdwLm6uuLixYuS\nZ7SwsEBhYaH29dWrVx/4M/eX2/3uPZaekpIClUoF4M4y7dixAzk5OdqvwsJCtG3bVjt+dSfiV6xY\ngfPnzyMuLg43b95ETEwMxJ098AdmbteuHVJTU7WvCwsLq92oVrZu7j/he29WV1dXrFmzRmfZbt26\nhT59+mD//v14//338dNPPyE3Nxc5OTmwsbGp1Ynh6vI86OIFlUqF5ORknfeSk5N1SpLkwVIwQl9/\n/TX27NmDli1bVhhWUFAAOzs7NGvWDHFxcfjhhx+0w6ZMmYJdu3bhp59+QllZGW7cuIG//voLQMNe\nZeLt7Y2oqCjk5OTg6tWr+Oijj+o9zaVLl6KoqAinTp3CunXrMHHiRADArFmz8N///ldbGtevX0dk\nZGSNp1tQUICWLVvCxsYG2dnZFU7MA1X/biZMmIBt27bh4MGDKCkpwaJFi6q9qqeydVPdxnfWrFlY\ntmwZTp8+DQC4efMmfvrpJwBAfn4+zMzM0KpVK5SUlOCtt95CXl5eldOqbBmqy9O6dWuYmJhU+UdE\ncHAwzp8/j4iICJSVlWHTpk04e/YsRo8eXWUGahwsBSPUsWNH+Pr6al/fu2H5/PPPsWjRIlhbW2PJ\nkiXajSdw5y/PqKgorFixAg4ODvDx8cGJEye007h/A1XdBqu6YVOnToWXlxfc3NwQFBSESZMmPfAv\nz3uH359FoVBg0KBBcHd3x9ChQ/Gf//xHe3XO/PnzERISguHDh8Pa2hp9+/ZFXFxcjXICwIIFC1BU\nVIRWrVrB398fwcHB1f4e7s3WrVs3fPbZZ5gyZQratWsHe3v7ag9VVbduKss6duxYvPLKK5g0aRJs\nbGzQo0cP7X0CQUFBCAoKQufOneHm5oaWLVtWOJxz/7Tvn351eZRKJV577TX069cP9vb2OHLkiM40\nHBwcsG3bNqxYsQKtWrXC8uXLsW3bNtjb21eZgRqHQjTkn3hERKTXJN1T0Gg08PHxwZgxYyodPm/e\nPHTq1AleXl6Ij4+XMgoREdWApKXw8ccfw9PTs9Jd8KioKFy4cAGJiYlYs2YNZs+eLWUUIiKqAclK\nIS0tDVFRUXjmmWcqPUkVGRmJ0NBQAICfnx9yc3N1LgMkIqLGZybVhJ9//nm8//77VV7RkJ6ernNS\nzdnZGWlpaTrXm/OZPEREdVPX08WS7Cls27YNjo6O8PHxqTbY/cMqK4G713zX5ausTIOHRr+L0jJN\nvaYj1debb74pewYuH5eNy6d/X0XFJegzdSXiz6ZXOrw+JCmFQ4cOITIyEh06dMDkyZOxZ88eTJs2\nTWcclUqlc+NOWlqa9oaihmJqagIbyxbIzS9q0OkSEcnpu+3x6Nm5Dby7tGvwaUtSCsuWLUNqaiou\nX76MjRs3YvDgwTpP0wSAkJAQ7XuxsbGwtbXVOXTUUOyslcjJYykQkWHIL7yNNT/H4sWpAyWZvmTn\nFO5197DQ6tWrAQBhYWEYOXIkoqKi4O7uDgsLC6xdu1aSedvbtER2XuGDR5RBQECA3BEkZcjLZ8jL\nBnD5mrKvf/0TA3w7onP71pJMv0nfvKZQKOp9fCxs6c8YN7g7gvy7NFAqIiJ5XM+5hRFzvsLWD0Ph\n0sa2yvHqs+00+Mdc2FkrkXOTh4+ISP998sMBjB/cvdpCqC+DLwV766Z7+IiIqKYupd1A1IGzeHai\nv6TzMfhS4IlmIjIE738bg5nj/WBnXfHpxg3J4Evhzp4CS4GI9Nex02k4kXgV00MelnxeBl8KdtYt\nkcPDR0Skp4QQeOebvXj+yQFo0dxc8vkZfCnY2yiRfZOlQET6aWdsIm4VlWBcYLdGmZ/hl4K1koeP\niEgvlZZp8N46NV55KgCmpo2zuTb4UrCzackTzUSkl36ITkDbVlYY9HDHRpunwZeCZctmKCktw+2S\nMrmjEBHV2M38YnwacRCvPTOkUZ8YbfCloFAoeFkqEemdTzcexLC+neDRwbFR52vwpQDcuSyVVyAR\nkb64nJ6NX/b8jReelOahd9UxilKw470KRKRH3vlmL2aO743WdhaNPm8jKQVelkpE+uHwiWScuXwN\nTz/aS5b5G0UptLazQFbuLbljEBFVS6Mpx5I1u/HqU4Fo3qxRPtmgAqMoBScHS2RmF8gdg4ioWpt3\nn4SlshlG9pfvUf/GUQr2lrh2g6VARE1XXkExVqzfh9dnNu4lqPczilJwtLfingIRNWkf/XAAQ/zc\n0bNTW1lzGEUpODlYIpN7CkTURJ1Luo6te0/hpWmD5I5iJKVgb4lr3FMgoiZICIHw1Tsxb3I/ONgo\n5Y4jTSkUFxfDz88P3t7e8PT0xMKFCyuMo1arYWNjAx8fH/j4+GDp0qVSRAEAWFk0R6lGg1tFJZLN\ng4ioLqIPnkNOXhGeGOUrdxQAgCTXPLVo0QJ79+6FUqlEWVkZ+vfvjwMHDqB///464w0aNAiRkZFS\nRNChUCi0ewsdVPaSz4+IqCaKikux7Ks9WPHiaJg10lNQH0SyFErlnd2gkpISaDQa2NtX3BgLIaSa\nfQWOPIRERE3Mqs2x8PFQwa+Hq9xRtCS7O6K8vBy+vr64ePEiZs+eDU9PT53hCoUChw4dgpeXF1Qq\nFZYvX15hHABYvHix9vuAgAAEBATUKY+jPe9VIKKmIzkjB99tO4btnz5d72mp1Wqo1er6hwKgEBL/\nuX7z5k2MGDEC7777rs4GPT8/H6amplAqlYiOjsb8+fNx/vx53XAKRYPtTSxZswttWllj5vjeDTI9\nIqK6EkLgqTd/Qt+ergh7rE+DT78+207JD2LZ2Nhg1KhROHr0qM77VlZW2kNMwcHBKC0tRXZ2tmQ5\nnByscC07X7LpExHVVPTBc7hyPQ9PyfR8o+pIUgpZWVnIzc0FABQVFWHnzp3w8fHRGSczM1PbZHFx\ncRBCVHreoaHcOafA5x8RkbwKCm9j6Ze7seTZ4Whmbip3nAokOaeQkZGB0NBQlJeXo7y8HFOnTsWQ\nIUOwevVqAEBYWBg2b96MVatWwczMDEqlEhs3bpQiihaff0RETcHHPxyAv5cb/Lo3nZPL95L8nEJ9\nNOQ5hUtpN/BM+Gbs+TKsQaZHRFRbZy5fw9TXNuL3Vc9IeqNakz6n0FTcvfqoCXcgERmw8nKB1z/7\nHS9OHdgk7lyuitGUgqWyORQKBQp4VzMRyeDHnSdQXi4wcYSX3FGqZTSlAPxzspkPxiOiRnY9pwDL\nv43BsrkjYGIi32Oxa8KoSsGJN7ARkQwWf7EL/xreEx4dneSO8kDGVQoOlsi8wXsViKjx7DqSiNOX\nMjFvcj+5o9SIUZUCn39ERI0pv/A23ly1E2/PDUKL5uZyx6kRoyoFJwcrXOU5BSJqJMvX70N/bzf4\ne7WXO0qNGVUptG1lhYzreXLHICIjcPxMOnYcPIeFMwLljlIrRlUKLk62SM3MlTsGERm426VlWPhp\nNN6YOQS2Vi3ljlMrRlUKzk42SM28KXcMIjJwKzcegmsbO4wa0FXuKLVmVKVga9UColwgr6BY7ihE\nZKD+vnAVEdEJeHvuCCgUTfuehMoYVSkoFAruLRCRZEpKNfjPR1H47zOD4WhvKXecOjGqUgAAZydb\npLEUiEgCn/94GO1aWWFcYDe5o9SZZB/H2VTd2VPgyWYialhnLmViw/bj2PbJU3p52Oguo9tTcHGy\n4Z4CETWo0jINXvpwO159OhBtWlnJHadejK4UnFkKRNTAVv14GI72lpgwpLvcUerN6A4fufDwERE1\noBOJGfhu+3H89rF+Hza6y+j2FFT/7Cnww3aIqL6Kb5fiheXbsOjfQ/X+sNFdRlcK1hYtYG5miuy8\nIrmjEJGee29dDDw7OmLMIE+5ozQYSUqhuLgYfn5+8Pb2hqenJxYuXFjpePPmzUOnTp3g5eWF+Ph4\nKaJUiucViKi+DiYkYcfBc3hrznC5ozQoSUqhRYsW2Lt3LxISEnDixAns3bsXBw4c0BknKioKFy5c\nQGJiItasWYPZs2dLEaVSd65A4nkFIqqbvIJivPxRFN6bH6x3zzZ6EMlONCuVdz6YuqSkBBqNBvb2\n9jrDIyMjERoaCgDw8/NDbm4uMjMz4eSk+8lEixcv1n4fEBCAgICAemdzacMb2Iio7hZ/sRNDertj\n4MMd5Y4CAFCr1VCr1Q0yLclKoby8HL6+vrh48SJmz54NT0/dY27p6elwcXHRvnZ2dkZaWlq1pdBQ\nVI42SEzJavDpEpHh26o+hYTzGdj2yXS5o2jd/wdzeHh4nacl2YlmExMTJCQkIC0tDfv27au0xe6/\nAqixLudycbJB6lUePiKi2km9mou3Vu/CJy+HQNmimdxxJCH51Uc2NjYYNWoUjh49qvO+SqVCamqq\n9nVaWhpUKpXUcQDcOXzEUiCi2ijTlGPB8t8w6/G+6O7eRu44kpGkFLKyspCbe2ejW1RUhJ07d8LH\nx0dnnJCQEKxfvx4AEBsbC1tb2wqHjqTi2tYW6dfzUFKqaZT5EZH++zTiIJQtzDFjbC+5o0hKknMK\nGRkZCA0NRXl5OcrLyzF16lQMGTIEq1evBgCEhYVh5MiRiIqKgru7OywsLLB27VopolSqubkZ2rW2\nRsrVHLi7tGq0+RKRfor7OxU/RCdg+6dPwcRE/+9aro5CNOFbexUKhWR3Hj8TvhmPD+uJEf6dJZk+\nERmGvIJijJz7DcLnDMeQ3u5yx6mR+mw7je6O5rs6OtvjYtoNuWMQURMmhMDLH0VhaJ9OelMI9WW0\npfCQiwNLgYiqtXbrUWRk5WHhjEC5ozQa4y0FZwdcTGUpEFHl4s+m4/MfD2Plq2PR3Nx4Hiht1KVw\nKS2bT0slogpy8orw3HuReGdeEFza2Modp1EZbSnYWbeEuZkJrufckjsKETUh5eUCL32wDcH9OmNY\nH+O7EMVoSwHgeQUiqujLX44gJ78IL08PkDuKLIy6FDryvAIR3ePPU6n46tc4fPrKozA3M5U7jiyM\nuhQecnbAxbRsuWMQURNw42Yh5v8vEv9bMAoqRxu548jGyEuB9yoQEaDRlOP55b9hbGA3BPZ6SO44\nsjLuUnBxwCUePiIyeh9u2I+SUg1emDpQ7iiyM+pScHa0QdbNQtwqKpE7ChHJZPv+s9iiPo3PFo6F\nmalRbxIBGHkpmJqaoJOrA84lXZc7ChHJ4Mzla1j0+e/44rVxcLBRyh2nSTDqUgAAzw5OOH0pU+4Y\nRNTIcvOLMGvpL3jj30MN+vMRaoul8JATTl+6JncMImpEZZpyPPfeVozw74yxgd3kjtOksBQ6OnJP\ngcjIvP9tDCBgtDeoVcd4nvJUBY8OjjifnIUyTTlPMhEZgciY04g+eBZbP5zO/+crYfS/EUtlczg5\nWOJyOm9iIzJ0py5mIvyLnVj9+gTYWbeUO06TZPSlAPAQEpExuJZdgLAlP+OtOSPg0cFR7jhNliSl\nkJqaisDAQHTr1g3du3fHJ598UmEctVoNGxsb+Pj4wMfHB0uXLpUiSo14duTJZiJDVlhcgplvbcbE\nIC+MGtBV7jhNmiTnFMzNzfHhhx/C29sbBQUFePjhhzFs2DB4eHjojDdo0CBERkZKEaFWPDs64Zst\nf8odg4gkUF4u8MLybejk2hpzJ/rLHafJk2RPoU2bNvD29gYAWFpawsPDA1euXKkwXlP5gJu7h4+a\nSh4iajjvrt2L3IJiLHsuCAqFQu44TZ5CSLwlTEpKwqBBg3Dq1ClYWlpq34+JicH48ePh7OwMlUqF\n5cuXw9PTUzecQoF/3/P64X++iIjo/x375+uuNaj7H92SlkJBQQECAgLw+uuvY+zYsTrD8vPzYWpq\nCqVSiejoaMyfPx/nz5/XDadQNNpf76FvbMLU0b4Y6tepUeZHRNLad+wSXvxgO356/0m4tbOTO06j\nqs+2U7Krj0pLSzFhwgQ8+eSTFQoBAKysrKBU3nnWSHBwMEpLS5GdLd9lod5d2iHhbMVDXESkf84l\nXccLK7bh8/+ONbpCqC9JSkEIgRkzZsDT0xMLFiyodJzMzP8/hh8XFwchBOzt7aWIUyO+HiocO5Mu\n2/yJqGFczynAjPDNeH3mEPTq5iJ3HL0jydVHBw8exIYNG9CzZ0/4+PgAAJYtW4aUlBQAQFhYGDZv\n3oxVq1bBzMwMSqUSGzdulCJKjfl0aYeTF67yzmYiPZZ3qxjTF/2IicN78plGdST5ieb6aMxzCgAw\nfPZX+ODF0XxiIpEeul1ShumLfkQn11YInz3MqK80apLnFPSRr4cKx8/yEBKRvtFoyrHg/UjY2yjx\nZthQoy6E+mIp3MO3qwrHeV6BSK8IIfDG538gv/A2PnhpNEx5+Lde+Nu7h69HO+4pEOmZDzbsx98X\nruKL18ajubnRP/i53lgK9+iockBewW1czymQOwoR1cC3kUexfd8ZfBP+OCyVzeWOYxBYCvcwMVHA\np2s7HD/D+xWImrrImNP4YvMRrF86Ca1sLeSOYzBYCvfx9VDh6OlUuWMQUTX2xF3AkjW7sPatx+Hs\nZCN3HIPCUriPf8/2OJiQLHcMIqrC/uOX8fJHUfhy0WPo6sbPRWhoLIX79OzcFmmZN3HjZqHcUYjo\nPodPJGPB+5FY9do4eHdpJ3ccg8RSuI+5mSl6d3fB4b+4t0DUlMT9nYq572zBZwvH8fEVEmIpVKKf\nd3sc/CtJ7hhE9I/4s+mYs+xXfPRyCPr0dJU7jkGrthRmzZoFIQTmzJnTWHmaBH8vNxxMSJI7BhEB\nOJl4FTPf+hnvPz8KA3w6yB3H4FVZCsnJyejfvz8effRR+Pv7IznZeA6ndG7fCkXFZUi9mit3FCKj\nduZSJp5e/BPemReEwF4PyR3HKFRZCmq1GklJSTh58iSSkpKgVqsbMZa8FArFP4eQjKcIiZqaUxcz\nEbroRyyeNQzD+nSWO47RqLIUQkNDkZycjNjYWKSkpCA0NLQxc8mun7cbDsRfljsGkVGKP5uO6Ys2\n4a05wzFqQFe54xiVah+dfeXKFbRr1w4ZGRlo27at9v3MzEw4OTlJH66RH519r8wb+Rgx52v8+f1z\nMDczlSUDkTE6cjIFz76zBe8/P4qHjOpIskdnt2t35zrgtm3bIicnB1999RWGDBkCb2/vOs1Mnzg5\nWKF9W1scPZ0mdxQio7H/+GXMWfYrPn45hIUgk2ofKVhYWIitW7ciIiICCQkJyMvLw5YtWzBgwIDG\nyieroX06YWdsIvr2bC93FCKDt+tIIl79OBpfvD6e9yHIqMo9hcmTJ6N79+6IiYnBggULcPnyZdjZ\n2SEgIACmpsZxOGWoXyfsOpIo2yEsImOxff9ZLPwkGl8vfpyFILMqS+HMmTNwdHSEh4cHPDw8jKYI\n7tXVrTXKy4FzydfljkJksDbvOonw1TuxfslEeHVu++AfIElVWQoJCQlYu3Ytbty4gcDAQAwYMAD5\n+fm4evXqAyeampqKwMBAdOvWDd27d8cnn3xS6Xjz5s1Dp06d4OXlhfj4+LovhUQUCgWG9emEXUcu\nyB2FyOAIIbDqp8P46PsDiHhnCjw6Sn/xCtWAqKE///xTvPDCC8LFxUX07du32nEzMjJEfHy8EEKI\n/Px80blzZ3H69GmdcbZv3y6Cg4OFEELExsYKPz+/CtOpRTzJHIi/LB5dsE7uGEQGRaMpF2+u+kME\nzflKXM3KkzuOwanPtrPGzz565JFHsGLFCiQlJeGdd96pdtw2bdpor1CytLSEh4cHrlzR/eCayMhI\n7b0Pfn5+yM3NRWZmZu0arRH07u6ClKu5SLt2U+4oRAbhdmkZ5v1vK84mXcem956Ak4OV3JHoHrX+\nQFMTExMMGjSoxuMnJSUhPj4efn5+Ou+np6fDxeX/Tyg5OzsjLS2twv0Pixcv1n4fEBCAgICA2kau\nF3MzUwT5d8G2mDOY9XifRp03kaHJu1WMsCW/wN6mJb59619o3oyfqdwQ1Gp1gz11QtI1UlBQgMce\newwff/wxLC0tKwwX913Vo1AoKoxzbynI5dFATyxetZOlQFQP17ILMH3Rj+jVzRmL/j0UpqZ8SHND\nuf8P5vDw8DpPq8q1cujQoXpdillaWooJEybgySefxNixYysMV6lUSE39/4+9TEtLg0qlqvP8pNTL\n0wV5t27jbNI1uaMQ6aULqVmY8NJ3GDWgKxbPGsZCaMKqXDPr16+Hr68vJk6ciHXr1tXoqqO7hBCY\nMWMGPD09sWDBgkrHCQkJwfr16wEAsbGxsLW1bZRHZ9SFiYkCYwZ5IFJ9Wu4oRHpn37FLmPTKD3j+\nif54dqJ/pUcEqOmo9tlHwJ37FaKjo/HHH38gNzcXgwcPRlBQEPr161flvQsHDhzAwIED0bNnT+0/\ngGXLliElJQUAEBYWBgCYO3cuduzYAQsLC6xduxa+vr664WR89tH9zlzKxMwlv2Df17NgYsJ/1EQ1\n8e1vx/DZpkP4bOFY3pTWiOqz7XxgKdyrsLAQe/fuRXR0NA4fPoxjx47VaaY11ZRKQQiBoGe/xltz\nhsOvOz/5iag6ZZpyvLV6F2JPpuCrRY/Bta2t3JGMSqOVQmNrSqUAAN9s+RMnL1zFhy+NkTsKUZOV\nV1CMZ9/dAlMTE3zySgisLVrIHcnoSPaUVNI1fkh37I67gOybhXJHIWqSkq7kYPyL38HdpRW+evMx\nFoIeYinUgq1VSwzr0wk/7/5b7ihETc6+Y5fw+H++w/SQh/Fm2FCY8QojvVTrtbZ//348++yzUmTR\nC1OCvBGpwyQNAAAVQElEQVQRHd+kDmsRyam8XODTjQfxn4+isHLhWDw5yvfBP0RNVo1uXjt+/Dgi\nIiLw448/okOHDpgwYYLUuZosXw8Vmpmb4fCJFPh78XMWyLjlFRTjhQ+2ISevCFs/DEWbVnxkhb6r\nshTOnTuHiIgIbNq0Ca1bt8bjjz8OIUSD3UqtrxQKBaaO9sU3W/9kKZBRO5t0DbPf/hUDH+6IzxeO\nQzNz43u8viGq8vCRh4cHjh8/jt9//x379u3Dc889Z5SfqVCZCUO6469zGbiQmiV3FCJZRMacxhML\nIzBvcj+EzxrGQjAgVZbCL7/8gpYtW2LgwIGYNWsWdu/ezePo/2jR3BzTRvviy5/j5I5C1Khul5Th\nzS924oPv9mHD25MwbnB3uSNRA3vgfQoFBQXaz2neu3cvpk2bhnHjxmH48OHSh2ti9yncKyevCIEz\nV+P3z2fw0b9kFC6m3sBz722FWzs7vPNcMGyseLlpU9VoN69lZ2dj8+bN2LhxI/bs2VOnGdZGUy4F\nAFj8xU60aG6GV58KlDsKkWSEEPhx5wn8b60aL00bhElBXnx+URPHO5plkn7tJkbPW4udX8xEK1sL\nueMQNbi8W8V4beXvSEzOwievhKBz+9ZyR6Ia4B3NMlE52uDRgG5Y9eNhuaMQNbiEc1cwet5a2Fq2\nwJYPp7EQjAT3FOrpek4Bhs36CtGfPY22razljkNUb6VlGnzxUyzWbzuGJc+OQJB/F7kjUS3x8JHM\n3l27F3kFt7HsuSC5oxDVS2JKFl5csQ221i3x3vxg/qGjp3j4SGZhE/rg90PncD75utxRiOpEoynH\nmp+PYNIr32NykDe+fetfLAQjxT2FBrIu8ih2H7mA9Usn8soM0itJV3Lw0gfbYGZmiv/NH8nPPjAA\n3FNoAp4Y6YPM7ALsjE2UOwpRjZSXC3z72zFMeHE9Rg3oih+WTWYhEPcUGtKB+CT8d2U0/vj8GbRo\nbi53HKIqnUu6jv+u3AEFgP8tGImOzg5yR6IGxBPNTciz72yBaxtbvPJUgNxRiCoovl2KTzcewsYd\nCXhx2kBMGuHNzxw3QE3y8NHTTz8NJycn9OjRo9LharUaNjY28PHxgY+PD5YuXSpVlEa1eNYwbN51\nAicTr8odhUjHgfgkBD37NZIzchD92QxMCfZhIVAFku0p7N+/H5aWlpg2bRpOnjxZYbharcYHH3yA\nyMjIqsPp4Z4CAPy652+s+SUOWz8M5dMjSXY3bhbi7S93I+5UGpbMGY7AXg/JHYkk1iT3FAYMGAA7\nO7tqx9HHDX5NjA3shnatrPBpxEG5o5ARKy3TYF3kUYyY/RVa2Vngj1UzWAj0QDX65DUpKBQKHDp0\nCF5eXlCpVFi+fDk8PT0rjLd48WLt9wEBAQgICGi8kHWkUCjw7vxgjJ63Dv283dCnp6vckcjI7I+/\njCVrdsPJ3hI/vDOZj6gwcGq1usE+AE3SE81JSUkYM2ZMpYeP8vPzYWpqCqVSiejoaMyfPx/nz5/X\nDaenh4/u2nfsEl79JBrbPnkK9jZKueOQEUjOyMHbX+3BuaTreH3mYAz168T7ZoxQkzx89CBWVlZQ\nKu9sKIODg1FaWors7Gy54khi4MMdMXqgB178YBs0mnK545ABKyi8jfe/jcG4F9bDp2s7/PHFMxjW\npzMLgWpNtlLIzMzUNllcXByEELC3t5crjmT+EzoIt0vKsOK7fXJHIQNUUqrBt5FHEThzDTKy8hG9\n8mnMfrwvmpvLdmSY9Jxk/3ImT56MmJgYZGVlwcXFBeHh4SgtLQUAhIWFYfPmzVi1ahXMzMygVCqx\nceNGqaLIytzMFCtfHYtHn/8WXTs4ImRQxfMmRLVVXi4QGXMaH2zYh4ecHbB+yb/g0dFJ7lhkAHjz\nWiM5cykTU1/fhFWvjUOvbi5yxyE9JYRAzLFL+N+6GDRvZoZXnwqAXw9eyEC6eEeznth//DJeWPEb\nvl/Gq0Go9o78nYKPNhzAtZwCvBwagOF9eRKZKsdS0CO/7vkby9fvw6b3noCzk43ccaiJE0Ig9mQK\nPv7hIDKu5+HZif4YP6Q7zEz5LEuqGktBz3z72zF89WscIt6dAmdHFgNVJITAgYQkfBpxEFk5tzB3\nkj9CArqxDKhGWAp6aF3kUXyz9U/8sGwK9xhISwgB9dFLWLnpEG7mF+O5yf4YPcADpiwDqgWWgp76\nNvIoVv98BGvD/4UubjzHYMxul5Zh697T+OrXOJiZmmDOv/oiuF8XlgHVCUtBj21Vn8LSL3fjs4Xj\n0Ls7r0oyNjl5RdgQdRzfbTsOz45OmDm+N/y92vMEMtULS0HP7T9+Gc8v/w2vPBWAx4f1lDsONYJL\n6dlYt/UoImNOY0TfzpgxrhevSKMGw1IwABdSszDzrZ8xuNdDePXpQJib8ZHbhqa0TIOdsYn4Pioe\n55OvY9IIb0wb44vWdpZyRyMDw1IwELn5RXhhxTbcLCjGJy+HQMUrkwxCRlYeInb8hU2//wW3tnZ4\nYpQPRvh35qMoSDIsBQNSXi7w1a9x+PKXI3gzbBhGD/SQOxLVQWmZBjHHLuHHP07gz1OpCAnohilB\n3ryggBoFS8EAJZy7gpc+2A6Pjo54M2woWtlayB2JHkAIgdOXruHn3SfxW8xpuLWzx4ShPTBmoAcs\nWjaTOx4ZEZaCgSq+XYoPNuzHz7tOYt7kfnhilC9vXmqCrmUXYMveU/hl998oKCrB+CHdMX5wd7i1\nq/6TB4mkwlIwcOeTryN89S5k3yzE4lnD+AC0JuB6TgGiD55D1IFzOH0pE0H+nTFhSA/06uYCExNe\nTkryYikYASEEog+ew9tf7YF3l7Z4bnI/dHVzlDuWUbmWfbcIzuJs0jUM6eWO4P5dMNC3I5o340lj\najpYCkaksLgEG7bH48tf4vCwhwpzJ/mju3sbuWMZJCEEzl6+hj1/XsTePy/iQuoNDO7tjpH9u2CA\nbwdePURNFkvBCBUVlyJiRwLW/HIE3To6Yfqjj6CflxsPXdRTYXEJDiUkY8+fF6E+ehHmZqYY3Nsd\ng3s9hN49XFgEpBdYCkbsdkkZNu86iQ1R8SgqLsGkEd54bFgPXq1UQ6VlGvx1PgOHTyQj9kQK/jqf\ngZ6d2iCw150i6Ohsz0dOkN5hKRCEEPjrfAZ+iI7H74fPw79ne4zs3xWBvR6CpbK53PGajJJSDU5f\nysSRkyk49Fcyjp9Jh2tbW/h7tUefnu3Ru7sLrPj7Ij3HUiAdeQXF+P3weUQfPIejp9LQp6crgvp1\nwQCfDmhtZzx7EEIIXLmeh4RzGUg4l474s1dw5vI1uLaxRe8ervDv2R5+PVxga9VS7qhEDarJlcLT\nTz+N7du3w9HRESdPnqx0nHnz5iE6OhpKpRLr1q2Dj49PxXAshXrLKyjGrrgL+P3QORw+kQJVa2v4\ne7VHP283PNLNGdYWLeSO2CA0mnIkZ+TgzOXrOJt0DWcvX8eJxAxoNOXw6aqCd5e28OmqQs9Obbjn\nRAavyZXC/v37YWlpiWnTplVaClFRUVi5ciWioqJw5MgRzJ8/H7GxsRXDsRQaVJmmHH9fuIoDCUk4\nGJ+Ev85noF1ra/Ts3BZendqgR6e2cHd1aNJFUVhcgpSruUi+kovkjBxcTL2Bs0nXkJhyA63slOjq\n5oiuHRzRpX1r9OzUBs5ONjwnQEanyZUCACQlJWHMmDGVlsKsWbMQGBiIiRMnAgC6du2KmJgYODk5\n6YZjKUiqtEyDxJQsnDifgb8SM/D3hau4lJYNS2VzuLs44CFnB7i2tUXbVtZo19oKbVpZwdHOUrIP\nfrldWobrObdwPbsA17JvITM7H9ezb+HqjXykZOQi+WoOcvOL4eJkg/bt7ODW1g4dVPbo6tYand1a\n81wA0T/qs+2U5fq69PR0uLj8/wfKODs7Iy0trUIpAMDixYu13wcEBCAgIKAREhoHczNTeHZ0gmdH\nJ0wK8gZw54F8GVl5uJSWjQupWUi5ehPHTqchIysfV67nISevCFYWzWFr2QK21i1ha9USFi2boUUz\nM7RoZobmzczQzNwUQgDlQkAIgfJygfJ//nu7tAwFhSUoKLyt/W9+YQnyC4tRVFyKVrYWcLS3RGs7\nSzja3/n+YQ8Vxg3ujvZtbdHGwYqfRkZ0H7VaDbVa3SDTku2i6/tbrKpd/HtLgaRnYqKAytEGKkcb\nDPDtUGF4maYcNwuKcTO/CLn5xcjJL0JhUQmKS8rufN0uQ0lpGUxMFDBRKKBQKLTfm5goYG5mCiuL\n5rBs2ezOf5XNYaVsDktlM1hbtOB9FkR1cP8fzOHh4XWeliyloFKpkJqaqn2dlpYGlUolRxSqJTNT\nEzjYKOFgo5Q7ChFJQJb98JCQEKxfvx4AEBsbC1tb20oPHRERUeOSZE9h8uTJiImJQVZWFlxcXBAe\nHo7S0lIAQFhYGEaOHImoqCi4u7vDwsICa9eulSIGERHVEm9eIyIyMPXZdvIyDiIi0mIpEBGRFkuB\niIi0WApERKTFUiAiIi2WAhERabEUiIhIi6VARERaLAUiItJiKRARkRZLgYiItFgKRESkxVIgIiIt\nlgIREWmxFIiISIulQEREWiwFIiLSYikQEZEWS4GIiLQkK4UdO3aga9eu6NSpE957770Kw9VqNWxs\nbODj4wMfHx8sXbpUqihERFRDZlJMVKPRYO7cudi1axdUKhV69eqFkJAQeHh46Iw3aNAgREZGShGB\niIjqQJI9hbi4OLi7u8PNzQ3m5uaYNGkStm7dWmE8IYQUsyciojqSZE8hPT0dLi4u2tfOzs44cuSI\nzjgKhQKHDh2Cl5cXVCoVli9fDk9PzwrTWrx4sfb7gIAABAQESBGZiEhvqdVqqNXqBpmWJKWgUCge\nOI6vry9SU1OhVCoRHR2NsWPH4vz58xXGu7cUiIioovv/YA4PD6/ztCQ5fKRSqZCamqp9nZqaCmdn\nZ51xrKysoFQqAQDBwcEoLS1Fdna2FHGIiKiGJCmFRx55BImJiUhKSkJJSQk2bdqEkJAQnXEyMzO1\n5xTi4uIghIC9vb0UcYiIqIYkOXxkZmaGlStXYsSIEdBoNJgxYwY8PDywevVqAEBYWBg2b96MVatW\nwczMDEqlEhs3bpQiChER1YJCNOFLgBQKBa9QIiKqpfpsO3lHMxERabEUiIhIi6VARERaLAUiItJi\nKRARkRZLgYiItFgKRESkxVIgIiItlgIREWmxFIiISIulQEREWiwFIiLSYikQEZEWS4GIiLRYCkRE\npMVSICIiLZYCERFpsRSIiEiLpSAjtVotdwRJGfLyGfKyAVw+YyZZKezYsQNdu3ZFp06d8N5771U6\nzrx589CpUyd4eXkhPj5eqihNlqH/wzTk5TPkZQO4fMZMklLQaDSYO3cuduzYgdOnTyMiIgJnzpzR\nGScqKgoXLlxAYmIi1qxZg9mzZ0sRhYiIakGSUoiLi4O7uzvc3Nxgbm6OSZMmYevWrTrjREZGIjQ0\nFADg5+eH3NxcZGZmShGHiIhqSkjgp59+Es8884z29XfffSfmzp2rM87o0aPFwYMHta+HDBkijh49\nqjMOAH7xi1/84lcdvurKDBJQKBQ1Gu/Odr/qn7t/OBERSUuSw0cqlQqpqana16mpqXB2dq52nLS0\nNKhUKiniEBFRDUlSCo888ggSExORlJSEkpISbNq0CSEhITrjhISEYP369QCA2NhY2NrawsnJSYo4\nRERUQ5IcPjIzM8PKlSsxYsQIaDQazJgxAx4eHli9ejUAICwsDCNHjkRUVBTc3d1hYWGBtWvXShGF\niIhqo85nIxpQdHS06NKli3B3dxfvvvtupeM899xzwt3dXfTs2VMcP368kRPWz4OWb+/evcLa2lp4\ne3sLb29vsWTJEhlS1s1TTz0lHB0dRffu3ascR5/X3YOWT5/XXUpKiggICBCenp6iW7du4uOPP650\nPH1dfzVZPn1ef0VFRaJ3797Cy8tLeHh4iFdffbXS8Wq7/mQvhbKyMvHQQw+Jy5cvi5KSEuHl5SVO\nnz6tM8727dtFcHCwEEKI2NhY4efnJ0fUOqnJ8u3du1eMGTNGpoT1s2/fPnH8+PEqN5r6vO6EePDy\n6fO6y8jIEPHx8UIIIfLz80Xnzp0N6v+9miyfPq8/IYS4deuWEEKI0tJS4efnJ/bv368zvC7rT/bH\nXBj6PQ01WT5Af6+0GjBgAOzs7Kocrs/rDnjw8gH6u+7atGkDb29vAIClpSU8PDxw5coVnXH0ef3V\nZPkA/V1/AKBUKgEAJSUl0Gg0sLe31xlel/Uneymkp6fDxcVF+9rZ2Rnp6ekPHCctLa3RMtZHTZZP\noVDg0KFD8PLywsiRI3H69OnGjikZfV53NWEo6y4pKQnx8fHw8/PTed9Q1l9Vy6fv66+8vBze3t5w\ncnJCYGAgPD09dYbXZf1JcqK5NhrqnoamqiY5fX19kZqaCqVSiejoaIwdOxbnz59vhHSNQ1/XXU0Y\nwrorKCjAY489ho8//hiWlpYVhuv7+qtu+fR9/ZmYmCAhIQE3b97EiBEjoFarERAQoDNObdef7HsK\nhn5PQ02Wz8rKSrsbGBwcjNLSUmRnZzdqTqno87qrCX1fd6WlpZgwYQKefPJJjB07tsJwfV9/D1o+\nfV9/d9nY2GDUqFE4evSozvt1WX+yl4Kh39NQk+XLzMzUtnlcXByEEBWODeorfV53NaHP604IgRkz\nZsDT0xMLFiyodBx9Xn81WT59Xn9ZWVnIzc0FABQVFWHnzp3w8fHRGacu60/2w0eGfk9DTZZv8+bN\nWLVqFczMzKBUKrFx40aZU9fc5MmTERMTg6ysLLi4uCA8PBylpaUA9H/dAQ9ePn1edwcPHsSGDRvQ\ns2dP7cZk2bJlSElJAaD/668my6fP6y8jIwOhoaEoLy9HeXk5pk6diiFDhtR726kQ+nzqnYiIGpTs\nh4+IiKjpYCkQEZEWS4GIiLRYCkREpMVSIKoFExMTTJ06Vfu6rKwMrVu3xpgxY2RMRdRwWApEtWBh\nYYFTp06huLgYALBz5044Ozvr3V2+RFVhKRDV0siRI7F9+3YAQEREBCZPnqzXD1UjuhdLgaiWJk6c\niI0bN+L27ds4efJkhYesEekzlgJRLfXo0QNJSUmIiIjAqFGj5I5D1KBkf8wFkT4KCQnBSy+9hJiY\nGFy/fl3uOEQNhqVAVAdPP/007Ozs0K1bN6jVarnjEDUYHj4iqoW7VxmpVCrMnTtX+x6vPiJDwQfi\nERGRFvcUiIhIi6VARERaLAUiItJiKRARkRZLgYiItFgKRESk9X8lyyilL4PD9gAAAABJRU5ErkJg\ngg==\n"
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    }
   ],
   "metadata": {}
  }
 ]
}